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This  report  summarizes  atomization  simulations  performed  to  investigate  the  fole 
of  this  process  in  liquid  rocket  engine  combustion  instabilities.  The  research  involves  the 
development  and  application  of  a  series  of  nonlinear  free-surface  models  based  on 
boundary  element  methods  (BEMs).  Models  have  been  developed  to  study  the  effects  of 
unsteady  injection,  transverse  acoustic  wave  interactions  (such  as  those  experienced  in 
the  F-1  engine),  longitudinal  acoustic  wave  interactions, ^and  the  dynamics  of  droplets  in 
the  presence  of  an  acoustic  wave.  Results  for  each  of  these  four  models  is  presented  in 
this  document. 


2>vici  CyC’AXjrFx’ 


xviiuiv  xj  8; 


RtAnflArd  form  .?QR  (Raw  aQ\ 


THIS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE 
COPY  FURNISHED  TO  DTIC 
CONTAINED  A  SIGNIFICANT 
NUMBER  OF  PAGES  WHICH  DO 
NOT  REPRODUCE  LEGIBLY. 


AFOSR  Contract  F49620-94-1-0151 


MODELING  OF  LIQUID  JET  ATOMIZATION  PROCESSES 


Stephen  D.  Heister,  Associate  Professor 

School  of  Aeronautics  and  Astronautics 

Purdue  University 

1282  Grissom  Hall 

West  Lafayette,  IN  47907 


1  August  1996 


Final  Technical  Report  for  Period  1  July  1994  -  30  June  1996 


Approved  for  Public  Release, 
Distribution  is  Unlimited 


Prepared  for: 

AFOSR/PKA 

Attn.  Dr.  Mitat  Birkin 

110  Duncan  Avenue,  Suite  B115 

Bolling  AFB,  DC  20332-0001 


Contents 

1  Summary 

2  Research  Objectives 

3  Status  of  the  Research 

3.1  Comparison  with  Experimental  Images  of  Jet  Profiles 

3.2  Effects  of  Unsteady  Injection  on  Atomization 

3.3  Simulation  of  F-1  Engine  Tangential-Mode  Instability 

3.4  Nonlinear  Droplet  Response  to  Acoustic  Excitation 

3.5  Coupled  Gas/Liquid  BEM  Simulations  . 

3.6  BEM  for  Viscous  Flows . 

4  Professional  Activities 
4.1  Technology  Transfer/Coupling  Activities 

5  References 

6  Appendix  A  -  Liquid  Jet  Visualization 

7  Appendix  B  -  Nonlinear  Jet  Modeling 

8  Appendix  C  -  Acoustic  Interactions  with  Liquid  Jet 

9  Appendix  D  -  Acoustic  Interaction  with  a  Droplet 

10  Appendix  E  -  BEMs  for  Two-Fluid  Flows 


'f  X 


1  Summary 

This  research  program  has  focused  on  enhancing  the  understanding  of  jet  atomization  processes  and  their 
contribution  to  combustion  instabilities  in  liquid  rocket  engines.  During  the  past  two  years,  progress  has  been 
made  in  understanding  the  role  of  the  atomization  process  in  tangential-mode  instabilities  in  the  F-1  engine. 
Boundary  element  simulations  developed  in  this  research  have  been  compared  with  experimental  results, 
showing  excellent  agreement-  Three  new  models  have  been  developed  during  the  effort.  The  first  model 
involves  a  fully-coupled  gas/liquid  simulation  of  an  injection  process  under  unsteady  chamber  conditions. 
The  second  model  involves  a  zonal  approach  to  a  fully  viscous  simulation,  wherein  boundary  layers  are 
treated  via  an  integral  formulation  and  boundary  elements  are  used  at  the  viscous/inviscid  interface.  The 
third  model  addresses  the  nonlinear  response  of  a  droplet  to  an  imposed  acoustic  perturbation  has  been 
studied.  Results  from  all  models  are  described  in  this  report. 

2  Research  Objectives 

The  understanding  of  the  complex  combustion  phenomena  present  in  liquid  rocket  engines  begins  with  the 
fundamental  process  of  fuel  and  oxidizer  jet  atomization.  Since  the  atomization  process  can  be  greatly  effected 
by  acoustic  disturbances^,  it  appears  as  a  primary  focus“  in  studies  involving  combusiton  stability.  For  this 
reason,  a  focused  research  effort  has  been  conducted  to  develop  models  capable  of  providing  quaniiiaiive 
information  regarding  atomization  processes  (in  both  steady  and  unsteady  chamber  environments). 

The  objective  of  this  research  has  been  to  develop  a  series  of  models,  incorporating  increasingly  complex 
physics,  to  assess  the  role  of  atomization  in  the  combustion  instability  process.  The  models  have  centered  on 
the  use  of  Boundary  Element  Methods  (BEMs)  as  a  means  to  provide  accurate  description  of  these  complex, 
nonlinear  processes  under  arbitrary  unsteady  conditions.  The  models  have  demonstrated  a  capability  to 
have  calculations  proceed  beyond  atomization  events.  While  the  basic  BEM  techniques  are  inviscid,  recent 
development  of  a  zonal  model  using  an  integral  method  for  boundary  layer  modeling,  permits  a  full  viscous 
capability.  This  model,  described  in  Section  3.6  of  this  report,  is  the  first  primary  atomization  model  to 
provide  accurate,  fully  nonlinear  treatment  of  atomization  processes  under  full-scale  Reynolds  numbers  con¬ 
sistent  with  actual  engine  conditions.  While  other  viscous  models  exist^’"^,  they  can  only  provide  simulations 
for  low  Reynolds  numbers  (of  the  order  of  a  few  hundred),  while  real  jets  are  typically  in  the  Re  ^  10^  —  10^ 
range. 

With  these  capabilities,  we  have  sought  to  address  the  sizes  of  droplets  produced  under  various  steady 
and  unsteady  injection  conditions.  In  addition,  the  response  of  droplets  to  acoustic  processes  is  another 
objective  of  this  work.  The  following  section  details  the  status  of  these  developments.  Section  4  provides  a 
description  of  professional  activities  (including  publications  and  student  theses)  associated  with  this  project. 

3  Status  of  the  Research 

Seven  major  tasks  have  been  accomplished  during  the  two  year  research  program.  Collaboration  with 
Dr.  Steven  Collicott  (also  in  the  School  of  Aeronautics  and  Astronautics)  has  enabled  the  generation  of 
high-quality  experimental  data  with  which  to  compare  atomization  models  (See  Section  3.1).  The  eflfect 
of  unsteady  injection  processes  on  the  atomization  of  a  liquid  jet  is  addressed  in  Section  3.2.  Section  3.3 
provides  a  description  of  the  role  of  the  atomization  process  in  tangential-mode  combustion  instabilities 
observed  during  the  development  of  the  F-1  engine.  In  Section  3.4,  the  response  of  a  droplet  to  an  acoustic 
perturbation  is  quantified,  while  Section  3.5  describes  the  fully-viscous  model  discussed  above.  Finally,  results 
from  a  fully-coupled  gas/liquid  model  are  described  in  Section  3.6  and  the  fully  viscous  model  highlighted 
in  Section  2  is  described  in  Section  3.7. 

3.1  Comparison  with  Experimental  Images  of  Jet  Profiles 

Under  the  direction  of  Professor’s  Heister  and  Collicott,  Mike  Moses  developed  an  experimental  setup  capable 
of  providing  high-resolution  imaging  of  low-speed  liquid  jets  during  the  atomization  process.  The  imaging 
system  magnifies  and  records  a  shadow  image  of  the  water  jet.  A  pair  of  achromatic  lenses  of  focal  length 
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Figure  1:  Jet  Profile  Comparison  Near  the  Pinch  Location,  (a)  Experiment  of  Moses,  (b)  BEM  Result, 
Wer=17.6,  Bo=:0.0109,  k  =  0.447. 


20  cm  and  diameter  4  cm  are  used  to  form  a  1:1  image  of  the  jet  in  the  intermediate  image  plane,  which  is 
then  imaged  onto  the  CCD  detector  at  the  desired  magnification.  A  variable-diameter  circular  iris  is  located 
between  the  achromats  and  controls  the  effective  F-number  of  the  imaging  system.  The  short-pulsed  spark 
source  assures  an  instantaneous  sample  of  the  jet  and  drop  shapes. 

The  edge  detection  scheme  used  in  this  project  is  an  improvement  of  the  method  of  Collicott,  et  al.^, 
which  was  applied  to  an  anamorphic  imaging  system  to  measure  growth  rates  of  waves  on  liquid  jets.  The 
additional  purpose  of  the  improved  method  is  to  investigate  droplet  sizes.  Therefore,  the  edge  detection 
procedure  is  applied  to  the  image  in  both  directions.  Programming  effort  is  required  to  sort  edges  in  droplets, 
and  to  assemble  the  data  from  horizontal  and  vertical  scans.  The  fundamental  operation  performed  in  these 
methods  is  to  convolve  the  image  greyscale  data  with  the  theoretical  edge-response  of  the  imaging  system. 
In  the  present  work,  the  combined  effects  of  diffraction  and  aberrations  are  approximated  by  assuming  a 
circular  Gaussian  impulse  response.  More  details  regarding  the  experimental  apparatus  can  be  found  in  the 
paper  attached  in  Appendix  A. 

A  simulation  using  our  finite-length  jet  code^  was  run  to  assess  the  accuracy  of  the  BEM  in  reproducing 
experimental  results.  Figure  1  provides  a  comparison  of  the  calculation  with  the  observation.  The  agreement 
between  the  BEM  simulation  and  the  experimental  measurements  are  excellent;  “recoil”  waves  generated  as 
a  result  of  the  pinching  process  are  resolved  within  the  calculation,  as  is  the  size  and  relative  spacing  of  main 
and  satellite  droplets.  This  comparison  confirms  the  high-accuracy  of  the  BEM  in  problems  of  this  nature. 

3,2  Effects  of  Unsteady  Injection  on  Atomization 

One  of  the  consequences  of  unsteadiness  in  the  combustion  chamber  is  the  possibility  of  inducing  unsteady 
massflow  through  injector  orifice  passages.  To  address  this  issue,  our  finite-length  jet  simulation^  was  applied 
to  a  range  of  unsteady  flows  in  which  the  frequency  and  amplitude  of  the  unsteady  component  of  injection 
velocity  was  varied.  The  effect  of  oscillation  magnitude  on  the  jet  breakup  is  investigated  by  six  simulations 
with  unsteady  velocity  component  {q')  varying  between  2  and  7%  of  the  mean.  Figure  2  shows  jet  profiles 
for  three  of  these  disturbance  magnitudes.  Each  pair  of  profiles  shows  the  jet  at  times  just  before  a  main 
droplet  and  a  satellite  droplet  are  shed  from  the  calculation.  As  this  figure  shows,  increasing  the  size  of  the 
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perturbation  decreases  the  breakup  length,  and  changes  the  shape  of  the  droplets. 

At  larger  disturbance  magnitudes,  the  main  droplets  take  on  a  '‘squashed”  shape  as  a  result  of  the  high 
velocity  fluid  encountering  lovi^er  velocity  fluid  which  has  already  exited  the  nozzle.  This  phenomena,  known 
as  the  "Klystron”  effect,  has  been  documented  qualitatively  by  numerous  researchers'’®  in  the  case  where 
forcing  perturbations  are  very  large  amplitude. 

Figure  3  shows  the  effect  of  dimensionless  frequency  {k)  on  the  character  of  the  jet  for  =  2%.  Each 
of  the  first  four  pairs  of  profiles  show  the  jet  at  times  just  before  a  main  droplet  and  just  before  a  satellite 
droplet  are  shed  from  the  calculation.  Increasing  the  frequency  from  k  =  0.5  to  /:  =  1.1  tends  to  decrease  the 
breakup  length  even  beyond  the  most  unstable  wave  number,  kmax  =  0.7;  a  trend  not  predicted  by  linear 
theory.  In  addition,  the  size  of  satellite  drops  tends  to  decrease  with  increasing  frequency  in  this  range.  At  a 
wave  number  of  ^  =  1.1,  satellite  drops  have  nearly  vanished  indicating  that  a  wave  number  near  this  value 
can  produce  monodisperse  atomization. 

Quantitative  predictions  of  droplet  sizes  are  shown  in  Fig.  4.  Here,  the  solid  lines  are  results  for 
infinitesimal  disturbances  and  experimental  results  are  a  composite  of  data  from  Lafrance^  and  Rutland  and 
Jameson^^.  The  finite-amplitude  results  are  for  g'  =  2%,  and  they  show  a  tendency  to  drive  conditions 
toward  a  monodisperse  result  by  increasing  satellite  drop  size  and  decreasing  main  drop  size  over  the  entire 
range  of  wave  numbers  investigated.  Once  again,  the  peak-sharpening  and  trough  broadening  which  occur 
with  increased  disturbance  magnitude  provide  an  explanation  for  the  observed  results.  At  higher  disturbance 
amplitudes,  local  curvature  in  the  transition  region  between  peaks  and  troughs  can  be  driven  to  large  enough 
values  to  promote  nonlinear  instability  at  wavelengths  predicted  to  be  stable  using  linear  analysis. 

A  complete  description  of  these  simulations  is  provided  in  the  Physics  of  Fluids  article  attached  in 
Appendix  B. 


(a) 


(b) 


(c) 


Figure  2:  Effect  of  Longitudinal  Disturbance  Amplitude  on  Behavior  of  Liquid  Jet  at  VFe  =  100,  Ug  =  uj  = 
0.7.  (a)  =  2%,  (b)  q'  =  4%,  (c)  g'  =  6% 
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(d) 


Figure  3:  Effect  of  Longitudinal  Disturbance  Frequency  on  Behavior  of  Liquid  Jet  at  We  =  100,  =  2%. 

(a)  u^/ujg  =  0.71,  (b)  uj/ujg  =  1.0,  (c)  uj/ufg  =  1.29,  (d)  uj/ujg  =  1.57 

3.3  Simulation  of  F-1  Engine  Tangential-Mode  Instability 

Through  the  recent  development  of  a  2-D  simulation  of  a  liquid  column  responding  to  an  imposed  crossflow, 
we  have  been  able  to  construct  simulations  pertinent  to  tangential-mode  acoustic  instabilities  in  liquid  rocket 
engines.  As  a  result,  efforts  have  been  focused  to  analyze  both  stable  and  unstable  injector  designs  from  the 
F-1  engine  test  program^ As  one  might  expect,  the  response  of  the  column  grows  dramatically  when  the 
acoustic  frequency  is  near  that  of  the  column  natural  frequency  (u;). 

Using  the  model,  a  simulation  was  conducted  of  a  highly-unstable  F-1  injector  configuration  (the  Double 
Row  Cluster,  DRC),  as  well  as  the  final,  stable  configuration  demonstrated  in  Flight  Rating  Tests  (FRT). 
The  DRC  design  exhibited  a  IT  instability  at  454  Hz  which  led  to  chamber  pressure  oscillations  of  the  order 
of  400%  of  the  mean.  The  most  prominent  difference  between  these  two  injectors  is  the  fuel  orifice  size  (3.57 
mm  radius  on  FRT  vs.  1.4  mm  radius  on  DRC).  The  combination  of  the  acoustic  frequency  in  the  chamber 
and  the  natural  frequency  of  the  DRC  fuel  colum  leads  to  conditions  very  near  the  resonant  frequency.  In 
fact,  we  calculate  that  ujg/ui  =  1.7  for  the  DRC,  while  the  FRT  design  has  ujg/u;  =  7.0,  a  value  far  from  the 
high  response  region. 

To  assess  the  impact  of  the  fuel  orifice  design  differences  between  the  two  injectors,  we  have  completed  a 
simulation  for  both  designs  at  a  fixed  Weber  number  of  0.1.  Results  of  the  column  shapes  at  various  times 
are  shown  in  Fig.  5  for  both  designs.  This  figure  shows  that  the  unstable  DRC  design  undergoes  violent 
oscillations,  whereas  the  stable  FRT  design  is  relatively  unaffected  by  the  imposed  oscillation.  Clearly,  the 
large  deformations  of  the  DRC  design  will  have  substantial  effect  on  the  jet  impingement  region  -  a  critical 
design  feature  for  this  impinging  element  injector.  For  this  reason,  we  believe  that  the  sensitivity  of 
the  DRC  design  to  transverse  acoustic  energy  may  be  a  major  contributor  to  the  instability 
of  this  injector  design. 

Further  description  of  this  analysis  can  be  found  in  the  Journal  of  Propulsion  and  Power  divticle  attached 
as  Appendix  C.  In  addition,  a  new  manuscript,  which  will  explore  this  problem  in  more  depth  is  currently 
in  preparation. 
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Figure  4:  Effect  of  Longitudinal  Disturbance  Frequency  on  Drop  Size  for  Liquid  Jet  at  We  =  100,  =  2%. 

Solid  Lines  are  for  q'  0. 

3.4  Nonlinear  Droplet  Response  to  Acoustic  Excitation 

A  model  was  developed  to  investigate  the  deformation  (and  possible  breakup)  of  a  droplet  exposed  to  an 
acoustic  wave.  A  series  of  simulations  at  a  constant  gas/liquid  density  ratio  of  0.001,  and  at  a  Weber  number 
{We)  of  0.58.  These  conditions  are  roughly  equivalent  to  a  100  micron  water  droplet  excited  by  a  160  decibel 
sound  wave.  Figure  6  shows  the  nonlinear  frequency  response  of  the  droplet  under  these  conditions.  Here, 
the  overall  level  of  response  is  inferred  from  the  aspect  ratio  (in  either  prolate  or  oblate  form)  of  the  droplet 
under  maximum  deformation  conditons.  The  droplet  response  is  much  more  complex  than  that  of  the  liquid 
column  in  that  4th-mode  coupling  is  present  in  many  cases.  Sharp  peaks  are  realized  at  several  higher-order 
harmonics.  The  region  0.9  <  <  1  is  characterized  by  actual  fragmentations  of  the  droplet. 

Figure  7  highlights  the  droplet  breakup  modes  identified  using  this  model.  At  lower  We  values,  small 
“nipples”  are  pinched  from  the  main  body  of  fluid.  As  We  grows,  dumbell-shaped  structures  are  encountered, 
followed  by  “doughnut-shaped”  structures.  At  the  very  high  We  values,  the  droplet  flattens  to  a  disk,  and 
small  rings  are  shed  from  the  periphery.  This  behavior  has  been  noted  by  several  experimentalists  by  exposing 
droplets  to  shock  waves,  thereby  causing  a  large  dynamic  pressure  about  the  drop.  Breakup  times  decrease 
dramatically  as  We  is  increased. 

A  complete  description  of  this  model,  which  includes  many  more  results,  is  provided  in  Appendix  D.  This 
manuscript  is  currently  in  review  in  International  Journal  for  Multiphase  Flows, 

3.5  Coupled  Gas/Liquid  BEM  Simulations 

During  the  peist  year,  we  have  developed  a  fully-coupled,  nonlinear  simulation  of  a  liquid  jet  issuing  into 
a  quiescent  gas.  This  model  has  been  developed  to  assess  unsteady  chamber  effects  on  the  atomization 
process.  Researchers  have  speculated  that  the  unsteady  chamber  conditions  existing  during  combustion 
instabilities  could  impact  the  atomization  process,  thereby  amplifying  (or  damping)  the  instability.  While 
several  calculations  have  been  made  using  our  model,  a  single  highlight  will  be  discussed  here.  Figure  8 
illustrates  the  influence  of  chamber  gas  density  on  the  behavior  of  the  jet  at  fixed  inflow  (i.e.  Weber  number) 
conditions.  The  gas  density  effect  is  measured  through  the  input  gas/liquid  density  ratio,  e.  Mushroom- 
capped  structures  appear  at  high  e  values  due  to  the  substantial  momentum  required  to  displace  the  “heavy” 
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Figure  5:  Simulation  of  Fuel-Jet  Response  in  the  Unstable  Double- Row  Cluster  (DRC)  and  Stable  Flight 
Rating  Tests  (FRT)  F-1  Engine  Injectors  (times  are  measured  in  seconds  from  start  of  disturbance) 
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Figure  7:  Droplet  Breakup  Modes  Encountered  Under  Acoustic  Excitation  at  the  Drop’s  Natural  Frequency 
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X 

Figure  8:  Effect  of  Gas  Density  on  Initial  Liquid  Jet  Behavior,  We  =  17.6 

gas  in  the  chamber.  Such  conditions  exist  in  high  pressure  combustion  devices  such  as  LREs. 

Presently,  we  are  investigating  the  injector  frequency  response  (both  amplitude  and  phase  shift)  under 
a  variety  of  unsteady  chamber  pressure  conditions.  Appendix  E  provides  some  information  in  this  regard, 
while  additional  results  will  be  presented  in  an  upcoming  manuscript. 

3.6  BEM  for  Viscous  Flows 

Substantial  progress  has  also  been  made  in  the  development  of  a  model  capable  of  addressing  liquid-ph2ise 
viscous  contributions  to  the  atomization  process.  Unfortunately,  we  were  forced  to  abandon  the  Dual  RecL 
procity  Method  for  which  we  presented  some  preliminary  results  last  year^^.  This  method  gave  fine  solutions 
for  steady  flows,  but  for  unsteady  flows  the  resulting  matrices  turn  out  to  be  ill-conditioned,  making  numer¬ 
ical  inversion  processes  inaccurate.  For  this  reason,  we  have  adopted  a  zonal  scheme  in  which  the  viscous 
region  is  treated  separately  via  an  integral  method  formulation.  Full  unsteady  simulations  are  possible  using 
this  methodology  within  the  framework  of  the  BEM  code  which  is  used  to  solve  for  properties  within  the 
inner,  inviscid  region. 

A  schematic  representation  of  the  zonal  approach  is  outlined  in  Fig.  9.  Our  original  BEM  is  used  to 
solve  for  conditions  (nodal  velocities  q  and  of  the  boundary  layer.  Using  4th-order  velocity  profiles  with 
appropriate  solid  wall  or  free  surface  conditions  at  the  outer  edge,  the  integral  method  provides  an  ordinary 
differential  equation  of  the  local  thickness,  S.  Continuity  and  conservation  of  momentum  (applied  at  the 
outer  boundary)  provide  relationships  for  the  surface  velocities,  u,  and 

Results  from  a  typical  calculation  are  shown  in  Fig.  10.  The  evolution  of  the  boundary  layer  and  free 
surface  is  depicted  for  a  case  where  fUe  —  10, 000,  Re  =  10, 000  and  the  initial  boundary  layer  thickness  (one 
radii  upstream  of  injection  point)  is  5%  of  the  orifice  radius.  These  are  typical  injection  conditions  for  actual 
rocket  injector  elements.  The  calculations  are  initiated  by  “ramping  up”  viscosity  to  the  desired  level  over 
a  period  of  2  time  units.  The  growth  of  the  boundary  layer  during  this  transient  generates  a  surface  wave 
which  convects  downstream.  At  long  times,  a  steady  solution  exists  even  under  these  high-speed  injection 
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Figure  9:  Schematic  Representation  of  Zonal  BEM/Integral  Method  for  Addressing  High-Speed,  Viscous  Jet 
Atomization  Processes 


conditions.  Obviously,  the  presence  of  the  gas  and  turbulence  effects  (which  are  not  presently  modeled)  can 
be  important  in  this  case.  Future  efforts  will  be  aimed  at  quantifying  interactions  due  to  these  physical 
processes. 


4  Professional  Activities 

The  efforts  outlined  in  the  previous  section  of  this  report  were  made  possible  by  two  grants  from  AFOSR,  A 
single  student,  Mr.  James  H.  Hilbing,  was  supported  under  the  base  grant  (F49620-94-1-0151).  In  addition, 
an  AASERT  grant  (F49620-93-1-0363)  was  utilized  to  support  Chris  A.  Spangler,  Mark  W.  Rutz,  Michael 
P.  Moses,  Ian  F.  Murray,  and  Kurt  Rump  (all  U.S.  citizens).  The  following  theses  were  written  as  a  result 
of  these  two  grants: 

Ph.D.  Dissertation 


Hilbing,  J.  H.,  ‘‘Nonlinear  Modeling  of  Atomization  Processes”,  August,  1996. 

M.S.  Theses 

Spangler,  C.  A.,  “Nonlinear  Modeling  of  Jet  Breakup  in  the  Wind-Induced  Regime”,  August,  1994. 
Rutz,  M.  W.,  “Effect  of  Transverse  Acoustic  Oscillation  on  the  Behavior  of  a  Liquid  Jet”,  December,  1995. 
Moses,  M.  P.,  “Visualization  of  Liquid  Jet  Breakup  and  Droplet  Formation”,  May,  1995. 

Murray,  I.  F.,  “Nonlinear  Modeling  of  the  Acoustically-Induced  Oscillations  of  Droplets”,  August,  1996. 
Rump,  K.  M.,  “Nonlinear  Model  for  a  Liquid  Jet  Injected  into  a  Quiescient  Gas”,  December,  1996. 
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A  list  of  journal  publications  (and  submissions)  associated  with  these  efforts  are  provided  in  the  following 
list.  Highlighted  items  (*)  have  been  attached  in  the  Appendices  of  this  report. 

Refereed  Journal  Publications  and  Submissions 


1.  Spangler,  C.  A.,  Hilbing,  J.  H.,  and  Heister,  S.  D..  “Nonlinear  Modeling  of  Jet  Atomization  in 
the  Wind-Induced  Regime”,  Physics  of  Fluids,  V  7,  No.  5,  pp  964-971,  1995. 

2.  Hilbing,  J.  H.,  Heister,  S.  D.,  and  Spangler,  C.  A.,  “A  Boundary  Element  Method  for  Atomization 
of  a  Finite  Liquid  Jet”,  Atomization  and  Sprays,  V  5,  No.  6,  pp  621-638,  1995. 

3.  *Hilbing,  J.H.,  and  Heister,,  S.D.,  “Droplet  Size  Control  in  Liquid  Jet  Breakup,”  Physics  of 
Fluids,  V  8,  No.  6,  pp.  1574-1581,  1996. 

4.  *Heister,  S.  D.,  Rutz,  M.,  and  Hilbing,  J.,  “Effect  of  Acoustic  Perturbations  on  Liquid  Jet  Atom¬ 
ization”,  To  Appear,  Journal  of  Propulsion  and  Power,  1995. 

5.  *Heister,  S.  D.,  “Boundary  Element  Methods  for  Two-Fluid  Free  Surface  Flows”,  In  Review, 
Engineering  Analysis  with  Boundary  Elements,  1996. 

6.  *Murray,  I.  F.,  and  Heister,  S.  D.,  “On  the  Response  of  a  Droplet  to  Acoustic  Excitation”,  In 
Review,  International  Journal  of  Multiphase  Flow,  1996. 

7.  Hilbing,  J.  H.,  and  Heister,  S.  D.,  “Coupled  Boundary  Element/Integral  Method  for  High-Speed 
Free-Surface  Flows”,  In  Preparation,  Engineering  Analysis  with  Boundary  Elements,  1996. 

8.  Rump,  K.  M.,  and  Heister,  S.  D.,  “Transient  Response  of  Liquid  Jet  under  Unsteady  Chamber 
Conditions”,  In  Preparation,  Atomization  and  Sprays,  1996. 

A  list  of  the  conference  papers  presented  in  association  with  work  under  these  grants  is  provided  in  the 
list  below.  The  starred  item  is  included  in  Appendix  A  of  this  report. 

Conference  Papers  and  Presentations 


1.  Hilbing,  J.H.,  and  Heister,,  S.D.,  “A  Boundary  Element  Method  for  Liquid  Jet  Atomization 
Processes”,  ILASS-94  Conference  Proceedings,  5  pages,  June  1994. 

2.  Spangler,  C.  A.,  and  Heister,  S.  D.,  “Nonlinear  Modeling  of  Jet  Atomization  in  the  Wind-Induced 
Regime”,  ILASS-94  Conference  Proceedings,  5  pages,  June,  1994. 

3.  Heister,  S.  D,,  Rutz,  M.,  and  Hilbing,  J.,  “Effect  of  Acoustic  Perturbations  on  Liquid  Jet  Atom¬ 
ization”,  AIAA  95-2425,  31st  AIAA  Joint  Propulsion  Conference,  San  Diego,  CA,  1995. 

4.  Hilbing,  J.H.,  and  Heister,,  S.D.,  “Developments  in  Nonlinear  Modeling  of  Atomization  Pro¬ 
cesses”,  ILASS-95  Conference  Proceedings,  1995. 

5.  *Moses,  M.  P.,  Collicott,  S.  H.,  and  Heister,  S.  D.,  “Visualization  of  Liquid  Jet  Breakup  and 
Droplet  Formation”,  7th  International  Symposium  on  Flow  Visualization,  Seattle,  WA, 

6.  Hilbing,  J.H.,  Heister,,  S.D.,  and  Rump,  K.,  “Recent  Advances  in  Nonlinear  Modeling  of  Atom¬ 
ization  Processes”,  ILASS-96  Conference  Proceedings,  1996. 

7.  Murray,  I.  F.,  and  Heister,,  S.D.,  “Modeling  Acoustically-Induced  Oscillations  of  Droplets”, 
ILASS-96  Conference  Proceedings,  1996. 

4.1  Technology  Transfer/Coupling  Activities 

Numerous  technology  transfers  have  occurred  during  the  period  associated  with  these  grants.  These  items  are 
summarized  the  the  table  provided  above.  Models  currently  under  development  should  be  of  great  interest 
to  the  liquid  rocket  engine  community,  since  we  soon  plan  to  have  a  fully  3-D  capability  for  high-speed  jets. 
Current  models  are  also  applicable  to  impinging  element  injectors  and  blanching  problems  associated  with 
oxidizer  deposition  on  chamber/injector  surfaces. 
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Table  1:  Atomization  Modeling:  Technology  Transfer  and  Transition 


Customer 

Result 

Application 

Lock  heed /Mart  in 

J.  Chrusciel 
(408)  756-3890 

Jet  Model  for  Projectile 

Control  Application 

Kigh-Speed  Interceptors 

Thiokol 

D.  Kawkins/T.  Boardman- 
(801)  863-3177 

Provided  Kybrid  Rocket  Ballistic 
Computer  Code 

HYBALID  Ballistics  Code 
Kybrid  Propulsion 
Demonstration  Program 

Lockheed/Martin 

Jim  Tegart 
(303)  977-9740 

Boundary  Element 

Modeling  of  Fluid 

Slosh  Dynamics 

Propellant  Tank 

Slosh  Analysis 

Kewlett  Packard 

Dr.  Graham  Ross 
(619)  487-4100 

Modeling  Liquid 

Jet  under  Excitation 

Ink  Jet  Printers 

Rockwell  Space  Systems 
Mark  Ventura 
(310)  922-0075 

Boundary  Element 

Modeling  of  Fluid 

Slosh  Dynamics 

Propellant  Tank 

Slosh  Analysis 
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Visualization  of  Liquid  Jet  Breakup  and  Droplet  Formation 

Michael  P.  Moses,  Steven  H.  Collicott,  and  Stephen  D.  Heister 
School  of  Aeronautics  and  Astronautics 
Purdue  University 
West  Lafayette,  IN  47907-1282 


Introduction 

Liquid  atomization  systems  are  commonplace  in  today’s  industry.  Some  of  the  many 
applications  include  fuel  injection  systems,  paint  application,  and  spray  drying.  Even  for 
relatively  low  speed  jets,  numerous  applications  exist.  Among  these  are  ink  jet  printers, 
droplet  generators,  and  manufacturing  processes  involving  the  formation  of  monodisperse 
droplets  from  a  column  of  molten  metal.  Due  to  the  large  number  of  practical  applications, 
the  literature  available  on  this  subject  is  quite  substantial.  Theoretical  treatments  of  this 
problem  date  back  to  the  late  1800’s  and  the  fundamental  linear  analysis  of  the  capillary 
instability  due  to  RayleighL  Linear  analyses  have  since  increased  in  sophistication  to  in¬ 
clude  the  presence  of  the  gas  phase^,  viscosity^,  and  three-dimensionality"*.  Unfortunately, 
these  models  can  not  address  important  physical  processes  which  occur  at  finite  surface 
deformations. 

Experiments®"®  have  revealed  that  a  given  wave  on  the  surface  of  the  jet  breaks  into  two 
droplets  rather  than  a  single  drop  as  postulated  by  the  linear  theory.  Rutland  and  Jameson®'^ 
examined  the  size  of  main  and  satellite  droplets  formed  during  this  nonlinear  process. 
These  data  were  collected  by  manually  measuring  the  size  of  the  drop  images  on  photographic 
plates.  Lafrance®  conducted  similar  experiments  on  drop  size,  using  electrostatic  separation 
of  main  and  satellite  droplets  and  subsequent  monitoring  of  fluid  volumes  collected.  While 
this  process  is  more  automated  than  that  used  by  Rutland  and  Jameson,  it  does  not  permit 
the  recombination  of  main  and  satellite  drops  (due  to  collisions  subsequent  to  the  pinch-off 
process)  which  may  occur^. 

Theoretical  efforts  in  recent  history  have  also  focused  on  the  nonlinearities  and  the  for¬ 
mation  of  main  and  satellite  droplets  from  a  single  instability  wave.  Weak  nonlinearities 
were  considered  by  including  higher-order  corrections  to  the  linear  theory®’*®.  Recently,  full 
nonlinear  simulations  of  the  entire  breakup  process  have  been  developed  through  the  use  of 
Boundary  Element  Methods  (BEMs).  Complete  nonlinear  simulations  have  been  developed 
for  an  infinite  inviscid  jet.  Models  have  been  developed  which  neglect**  and  include*^  the 
presence  of  the  gas  phase.  Hilbing,  et  al.*®  included  the  orifice  geometry;  a  simulation  for  a 
finite  liquid  jet. 

At  present,  additional  measurements  of  droplet  sizes  are  required  for  comparison  with 
the  fully- nonlinear  models.  In  the  past,  measurements  of  droplet  sizes  from  images  has  been 
performed  manually,  limiting  the  data  available  on  drop  sizes.  With  digital  image  processing 
applied  to  the  data  collection  process,  a  large  amount  of  drop  size  data  can  be  processed 
efficiently  and  accurately.  Additionally,  measurements  on  the  influence  of  the  gas-phase 
pressure,  in  the  “first  wind-induced  regime”  of  break-up  are  sparse.  Previous  observations*’’® 
did  not  attempt  to  address  this  effect,  which  is  difficult  to  discern  at  low  jet  velocities.  Only 
with  accurate  resolution  of  droplet  sizes  and  wave  shapes  does  the  potential  exist  to  assess 
aerodynamic  effects  in  the  low-speed  region.  For  these  reasons,  research  was  initiated  to 
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expand  the  drop  size  database  and  to  address  aerodynamic  effects  within  the  low  speed 
regime. 

Experimental  Apparatus 

The  experimental  apparatus  is  a  water  tank  with  acoustic  perturber,  an  orifice,  and  the 
optical  system.  The  water  tank  is  a  low  pressure  thin- walled  steel  tank  56  cm  long  with 
a  12.5  cm  inner  diameter  and  spherical  end  caps.  The  head  height  of  water  is  fixed  to  be 
30.2  cm.  The  tank  is  mounted  on  a  vertical  traverse  with  a  range  of  approximately  115  cm, 
allowing  the  optical  setup  to  remain  at  a  fixed  height  on  the  table.  The  traverse  is  mounted 
to  an  optics  table  with  vibration-isolation  legs,  isolating  the  experiment  from  vibrations  in 
the  floor.  This  is  found  to  be  important  to  the  success  of  the  experiment.  The  low  pressure 
tank  is  used  with  an  orifice  of  diameter  0.0.566cm  and  L/D  =  5.11.  The  upper  chamber  of 
this  orifice  produces  a  contraction  ratio  of  approximately  120:1.  Distilled  water  is  used  for  all 
experiments,  with  water  temperature  measured  during  each  run  (fairly  constant  near  20°  C 
throughout).  Perturbations  are  introduced  into  the  liquid  jet  with  a  piezoelectric  speaker 
mounted  inside  the  top  of  the  tank. 

The  optical  system  is  shown  in  Fig.  1  and  consists  of  the  illumination  and  imaging  sys¬ 
tems.  Collimated  white  light  from  a  short-pulsed  (20  nsec)  spark  source  is  used  to  illuminate 
the  water  jet.  The  center  post  of  the  spark  source  casts  a  shadow  in  the  center  of  the 
beam.  Therefore,  the  centerline  of  the  imaging  system  is  offset  from  the  centerline  of  the 
illumination  system,  as  indicated  in  Fig.  1. 

The  imaging  system  magnifies  and  records  a  shadow  image  of  the  water  jet.  Because 
of  the  magnification  required,  the  object-to-camera  distance  is  small,  placing  the  lens  too 
close  to  the  jet  and  subjecting  it  to  wetting.  Therefore,  a  pair  of  achromatic  lenses  of  focal 
length  20  cm  and  diameter  4  cm  are  used  to  form  a  1:1  image  of  the  jet  in  the  intermediate 
image  plane,  which  is  then  imaged  onto  the  CCD  detector  at  the  desired  magnification.  A 
variable-diameter  circular  iris  is  located  between  the  achromats  and  controls  the  effective 
F-number  of  the  imaging  system. 

A  Super-VHS  format  CCD  camera  is  used  with  the  35  mm  lens  and  extension  tubes  to 
image  the  jet  break-up  process.  The  camera  has  a  0.64x0.48  cm  chip  with  a  resolution  of 
768  streamwisex494  transverse  pixels.  The  short-pulsed  spark  source  assures  an  instanta¬ 
neous  sample  of  the  jet  and  drop  shapes.  A  typical  image  is  shown  in  Fig.  2a. 

The  magnifications  of  the  complete  imaging  system,  i.e.,  from  object  space  to  pixels  in 
computer  memory,  are  determined  for  each  test  run  by  acquiring  images  of  a  target  ruler 
marked  in  0.5  mm  increments.  Magnifications  parallel  with  and  perpendicular  to  the  jet  axis 
must  be  measured  separately  for  the  CCD  camera  in  use.  Magnifications  are  determined  from 
spectral  analysis  (1-D  power  spectra)  of  numerous  lines  within  these  calibration  images.  The 
standard  deviation  of  the  magnification  measurements  is  found  to  be  approximately  0.5%  in 
all  cases.  Without  the  spectral  method,  magnification  measurement  was  found  to  give  3-5% 
error,  and  was  the  largest  error  source  in  the  experiment. 

Digital  Image  Processing 

The  edge  detection  scheme  used  in  this  project  is  an  improvement  of  the  method  of 
Collicott,  et  al.'^*,  which  was  applied  to  an  anamorphic  imaging  system  to  measure  growth 
rates  of  waves  on  liquid  jets.  The  additional  purpose  of  the  improved  method  is  to  investigate 
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droplet  sizes.  Therefore,  the  edge  detection  procedure  is  applied  to  the  image  in  both 
directions.  Programming  effort  is  required  to  sort  edges  in  droplets,  and  to  assemble  the 
data  from  horizontal  and  vertical  scans. 

The  fundamental  operation  performed  in  these  methods  is  to  correlate  the  image  greyscale 
data  with  the  theoretical  edge-response  of  the  imaging  system.  In  the  present  work,  the 
combined  effects  of  diffraction  and  aberrations  are  approximated  by  assuming  a  circular 
Gaussian  impulse  response.  This  impulse  response  is  analytically  integrated  in  one  direction 
to  determine  the  edge  response.  The  width  of  the  Gaussian  is  found  empirically  in  the  present 
work.  With  complete  specifications  of  the  lenses,  the  width  could  be  computed.  Empiricism 
is  found  to  suffice  for  the  present  work,  in  that  the  basic  functions  of  the  correlation:  peak 
formation,  noise  reduction,  and  process  automation,  are  not  highly  sensitive  to  the  choice 
of  the  width.  For  the  present  work,  the  speed  of  direct  computation  in  the  spatial  domain 
is  sufficient,  and  is  simpler  to  implement  in  sub-regions  of  an  image  than  an  FFT-based 
computation.  Each  1-D  correlation  is  scanned  for  peaks  caused  by  sharp  greyscale  ramps 
occurring  at  edges.  The  upper  plot  in  Fig.  3  shows  the  raw  intensity  levels  along  the  scan 
line.  The  black  portions  (/(x)  =  0)  are  the  drops  while  the  white  area  (/(x)  ~  250)  is  the 
background.  The  lower  plot  in  Fig.  3  is  the  result  of  the  correlation  of  the  scan  line  with  the 
edge  response. 

Only  peaks  with  intensity  greater  than  a  threshold  value  are  considered  valid  edge  loca¬ 
tions.  In  this  manner,  gradual  greyscale  shifts  (such  as  the  dark  region  at  the  top  of  Fig.  2) 
are  discarded.  Edge  detection  at  pixel  resolution  is  obtained  by  noting  the  pixel  location  of 
the  correlation  peak.  However,  edges  don’t  necessarily  lie  at  integral  pixel  locations,  so  a 
“center  of  mass”  calculation  yields  a  sub-pixel  edge  location.  The  correlation  method  reduces 
the  effects  of  image  noise,  as  seen  in  Fig.  3  where  the  correlation  is  smoother  than  the  image 
intensity. 

A  trapezoidal  integration  is  used  to  compute  the  area  of  the  droplet  image  assuming 
that  the  calculated  area  represents  the  maximum  cross-sectional  area  of  an  axisymmetric 
object.  Note  that  multiple  views  would  be  required  to  resolve  3-D  effects;  assumption  of 
axisymmetry  is  deemed  to  be  adequate  for  this  study. 

Data  Reduction  Validation 

Three  validation  tests  assess  the  accuracy  of  the  image  processing  and  drop  sizing  techniques. 
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Figure  2:  (a)  Greyscale  Image  of  Typical  Droplet  Formation,  (b)  BEM  Simulation  of  Ex¬ 
periment. 


The  first  work  compared  measured  and  theoretical  wave  growth  rates  in  the  Rayleigh  break¬ 
up  regime  ,  and  showed  that  the  concept  of  correlating  the  image  with  the  theoretical 
edge  response  is  an  accurate  method  of  edge  location.  That  demonstration  was  actually  an 
anamorphic  imaging  system,  with  a  streamwise  magnification  equal  to  1/40  of  the  transverse 
magnification. 

The  error  introduced  into  the  digital  image  processing  by  the  discretization  of  the  edge 
of  a  circle  into  N  points  was  also  considered  as  a  validation  exercise.  Results  indicate  that 
this  error  is  rapidly  reduced  as  the  number  of  points  is  increased,  and  that  for  N  >  75,  the 
error  is  less  than  0.1%.  Therefore,  to  eliminate  numerical  integration  error,  images  need  to 
be  magnified  such  that  the  small  satellite  drops  will  have  at  least  75  pixels  on  the  perimeter. 

The  third  verification  test  examines  the  accuracy  of  the  entire  data  acquisition  and  re¬ 
duction  process  in  determining  the  size  of  droplets.  A  drill  gauge  plate  is  placed  in  the  object 
plane,  and  images  of  the  holes  are  acquired.  Error  analysis  shows  that  there  is  a  ±0.7%  error 
in  drop  size. 

Results 

Experiments  are  performed  at  two  Weber  numbers  in  order  to  keep  jet  velocity,  and  hence, 
the  aerodynamic  effects,  fixed  for  different  perturbing  frequencies.  The  Weber  number. 
We  =  pU^a/a,  is  the  dimensionless  parameter  characterizing  this  process.  Here,  U  is  the 
mean  orifice  exit  velocity  (from  volume  flowrate  measurements),  a  is  the  orifice  radius,  and 
p  and  <7  are  liquid  density  and  surface  tension,  respectively.  All  experiments  are  conducted 
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Figure  3:  Intensity  and  Correlation  Plots  for  Sample  Droplet  Image 

with  water  exhausting  into  ambient  pressure  air.  The  two  Weber  numbers  are  17.6  and  76.6, 
corresponding  to  jet  velocities  of  2.13  and  4.44  m/s,  respectively,  for  this  orifice. 

By  varying  the  frequency  of  the  speaker  at  the  top  of  the  tank,  various  dimensionless 
wave  numbers,  k,  are  generated.  Linear  stability  theories  indicate  that  the  jet  should  be 
unstable  for  0  <  A;  <  1,  but  in  practice  it  is  difficult  to  force  instability  for  very  long  (small 
k)  or  very  short  waves  due  to  the  tendency  of  the  jet  to  break  up  at  the  most  destructive 
wavenumber,  k  «  0.7.  Thus,  experimental  observations  are  limited  to  0.3  <  /t  <  0.85, 
similar  to  that  achieved  by  others^’®.  Frequencies  in  the  range  of  300-1800  Hz  produce  this 
wave  number  range  for  the  two  Weber  numbers  noted  above.  The  measured  wave  numbers, 
based  on  a  calculation  of  the  volume  of  fluid  within  a  single  wave  near  the  breakup  point, 
differ  noticeably  from  the  expected  wave  numbers.  Measured  wave  numbers  were  as  much  as 
10%  smaller  than  the  input  wave  number,  attributed  to  jet  (and  wave)  stretching  by  gravity 
and  surface  tension  effects. 

A  typical  image,  ai  k  =■  0.447,  and  We  =  17.6,  is  compared  with  numerical  simulation^® 
in  Fig.  2.  The  image  shows  the  main  and  satellite  droplets,  plus  a  wavy  surface  on  the 
core  column  of  liquid.  This  wavy  surface  is  another  satellite  droplet  just  prior  to  pinch. 
The  waves  are  caused  by  the  “recoil”  of  the  surface  after  the  main  drop  separates.  Note 
that  these  features  are  predicted  correctly  by  the  inviscid  boundary  element  method  (BEM) 
simulation,  as  are  the  sizes  of  the  main  and  satellite  droplets.  The  BEM  simulation  predicts 
atomization  of  the  satellite  as  well  into  two  subsatellites  whose  trajectories  merge  (indicating 
possible  coalescence).  Viscosity  may  serve  to  inhibit  this  process,  which  is  not  observed  in 
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this  particular  experiment  (but  can  be  observed  at  lower  wave  numbers  where  satellite  drops 
are  larger). 

A  summary  of  the  droplet  size  measurements  from  both  the  present  and  previous^* 
studies  is  compared  with  BEM  predictions^^  in  Fig.  4.  Wave  numbers  are  adjusted  such 
that  the  total  volume  of  fluid  in  main  and  satellite  drops  corresponds  to  the  initial  column 
(of  radius  a)  of  the  same  volume.  Each  of  the  data  from  this  study  represents  the  average 
of  between  two  and  thirteen  independent  measurements.  At  a  given  set  of  conditions,  the 
standard  deviation  of  main  drop  sizes  is  typically  less  than  2%,  while  satellite  drops  have  a 
larger  dispersion  (<  9%)  due  to  their  lesser  size. 

Agreement  of  the  present  experimental  results  with  those  of  Rutland  and  Jameson  is  quite 
good,  with  the  exception  of  four  satellite  drops  in  the  0.3  <  k  <  0.4  range.  Discrepancies  here 
could  be  caused  by  formation  of  subsatellites  (which  combine  with  main  drops).  Agreement 
of  the  present  data  with  measurements  of  Lafrance  indicates  a  systematic  bias  which  may  be 
attributable  to  the  electrostatic  data  collection  of  Lafrance.  In  that  technique,  one  cannot 
account  for  the  wave  stretching  phenomena  observed  in  this  study  and  by  Rutland  and 
Jameson. 

Agreement  of  the  present  results  with  the  theoretical  BEM  predictions  also  indicates  a 
systematic  bias.  This  is  believed  to  the  effect  of  the  inviscid  assumption  in  the  BEM  model. 
The  BEM  results  do  predict  an  increase  in  the  size  of  satellite  drops,  albeit  minor,  with 
increasing  Weber  number.  Comparison  of  present  experimental  results  indicates  a  similar 
trend,  but  the  effect  appears  to  be  more  pronounced  than  that  predicted  with  the  inviscid 
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Figure  5:  Jet  Surface  Shape  at  Different  Downstream  Positions,  We  =  17.6,  k  =  0.808 


theory.  Nevertheless,  the  overall  comparison  between  the  current  measurements  and  the 
inviscid  theory  is  quite  good. 

Numerous  images  of  wave  shapes  in  various  stages  of  the  pinch-off  process  are  also  ob¬ 
tained  in  this  work.  A  sampling  of  these  waves  is  shown  in  Fig.  5,  where  curves  (a)-(c) 
were  obtained  at  successive  downstream  locations.  Note  the  nonlinear  wave  growth,  and  the 
preferential  necking  of  the  surface  toward  the  downstream  (right)  half  of  a  given  wave.  The 
preferential  necking  can  not  be  predicted  with  an  infinite  jet  analysis,  so  simulations  of  a 
finite-length  jet  are  underway.  Many  other  results  of  jet  surface  and  droplet  images  in  this 
experiment  are  published  in  Moses^®. 

Conclusions 

Experiments  measuring  main  and  satellite  drop  sizes  show  the  influence  of  aerodynamic 
forces  as  Weber  number  is  increased,  causing  the  satellite  drops  to  become  larger.  This 
effect  is  also  predicted  by  the  BEM  modeP^.  This  effect  can  not  be  seen  in  the  previously 
available  data,  due  to  the  large  scatter  in  Weber  number®’®.  Investigation  of  wave  shape 
profiles  at  We=17.6  and  76.6  compare  very  favorably  with  BEM  models.  This  research  (see 
also  Moses^®)  has  also  expanded  the  database  of  experimental  wave  shape  profiles  available 
for  comparisons. 

The  benefits  of  the  1-D  theoretical  edge-response  correlation  are  sizable:  unambiguous 
and  consistent  edge  location,  a  large  reduction  in  the  effect  of  image  noise  on  edge  location, 
automation,  and  the  ease  of  error  quantification.  Wave  shape  determination  is  simple  and 
droplet  sizing  requires  only  additional  processing.  Magnification  measurement  by  spectral 
analysis  of  test  images  reduces  the  magnification  error  to  a  level  comparable  with  other  errors 
in  the  system.  These  improvements  make  possible  the  detailed  investigations  of  drop  sizes 
and  wave  shapes  for  a  wide  range  of  jet  flows. 
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A  Boundary  Element  Method  (BEM)  has  been  developed  to  investigate  the  nonlinear  evolution  of 
the  surface  of  liquid  Jets  injected  from  circular  orifices  under  unsteady  inflow  conditions.  For  fixed 
wavelength  perturbations,  the  model  predicts  the  formation  of  main  and  satellite  droplets.  The  size 
of  the  droplets  is  affected  by  changes  in  the  perturbation  wavelength,  perturbation  magnitude  and 
Weber  number.  Satellite  droplet  velocities  are  less  than  main  droplet  velocities  due  to  the  sequential 
shedding  of  droplets  from  the  orifice.  Using  this  information,  one  can  determine  the  likelihood  of 
droplet  recombinations  downstream  of  the  initial  pinching  event.  ©  1996  American  Institute  of 
Physics.  [S 1 070-663 1  (96)02506-8] 


I.  INTRODUCTION 

The  breakup  of  a  low  speed  liquid  jet  is  one  of  the  fun¬ 
damental  problems  in  two-phase  flow.  The  resultant  drop 
sizes  formed  from  the  breakup  process  are  of  interest  to  a 
variety  of  industrial  processes,  such  as  the  manufacture  of 
metal  powders  via  solidification  of  a  stream  of  molten  mate¬ 
rial.  In  addition,  the  performance  of  ink  jet  printers  is  closely 
tied  to  the  atomization  process  in  the  stream  of  ink  emanat¬ 
ing  from  an  injection  orifice.  For  these  reasons,  a  wide  vari¬ 
ety  of  both  analytic  and  experimental  treatments  have  been 
applied  to  this  problem. 

Early  analytic  works  were  performed  on  linearized  gov¬ 
erning  equations;  results  led  to  a  prediction  that  a  single 
droplet  is  formed  from  a  given  wavelength  instability  on  the 
surface  of  the  jet.  Experimental  observations,  and  more  re¬ 
cent  nonlinear  analytic  treatments,  have  revealed  that  a  given 
wave  along  the  surface  actually  atomizes  into  two  droplets 
(termed  “main”  and  “satellite”  drops  in  the  literature). 
Many  experiments  have  been  conducted  utilizing  a  small 
amplitude  perturbation  in  order  to  create  disturbances  of  a 
given  wavelength'”*  and  to  measure  the  size  of  resulting 
main  and  satellite  droplets.  Nonlinear  models  have  been  de¬ 
veloped  to  analyze  this  situation  by  assuming  a  periodic  dis¬ 
turbance  along  the  length  of  the  jet  (infinite  jet 
assumption),^  '^  In  addition,  the  effect  of  the  presence  of  the 
orifice  has  been  included  in  recent  “finite  length”  jet 
simulations.^  These  models  have  been  quite  successful  in 
replicating  the  droplet  sizes  obtained  from  the  experiments 
using  a  small  controlled  disturbance. 

In  spite  of  these  advances,  many  applications  would  ben¬ 
efit  from  an  ability  to  provide  additional  droplet  size  control 
via  modulation  of  disturbance  amplitude  or  frequency  from  a 
controlled  perturbation.  For  example,  wavelengths  which 
would  be  predicted  to  be  stable  on  a  linear  basis  could  be 
forced  to  instability  through  the  use  of  high  amplitude,  non¬ 
linear  perturbations.  Monodisperse  droplet  trains  have  been 
produced  in  this  fashion  by  Dressier.*'^  Another  alternative  is 
to  utilize  finite  amplitude  frequency  modulation  to  affect 
nonlinear  wave  deformation  processes  and  control  drop  size. 
These  notions  have  been  utilized  by  Orme  and  coworkers  in 
recent  experiments. This  technique  relies  on  recombina¬ 
tion  of  main  and  satellite  drops  due  to  “microspeed  disper¬ 


sions  between  the  two  bodies.  Through  experimentation,  ’ 
various  manufacturing  operations  have  made  use  of  these 
ideas  to  optimize  droplet  production  processes. 

While  experimentalists  and  manufacturing  disciplines 
have  made  progress  in  controlling  droplet  sizes  using  the 
techniques  described  above,  the  approaches  have  received 
very  little  analytic  treatment.  Models  based  on  the  infinite  jet  J 
treatment  cannot  be  used  to  analyze  these  disturbances  be-  I 
cause  of  the  assumption  of  periodicity,  i.e.  the  effect  of  dis-  1 
turbance  magnitude  increasing  as  one  approaches  the  orifice  ' 
cannot  be  addressed. 

However,  with  the  use  of  a  finite-length  jet  model.’  the 
effect  of  large  amplitude  perturbations  on  the  behavior  of  the 
jet  (and  resulting  drop  sizes)  can  be  predicted.  In  this  paper, 
we  discuss  the  model  developed  for  this  purpose  and  present 
results  for  a  variety  of  different  injection  conditions.  The 
model  is  described  in  Section  II.  Validation  of  the  model  is 
presented  in  Section  III,  while  results  are  described  in  Sec¬ 
tion  IV. 


II.  MODEL  DEVELOPMENT 


We  assume  an  axisymmetric,  incompressible,  inviscid 
flow  in  which  gas  pressure  variations  are  negligible.  Under 
these  assumptions,  the  unsteady  liquid  flow  is  described  by 
Laplace’s  equation  V’d)=0,  where  (f>  is  the  velocity  poten¬ 
tial  defined  as  the  function  whose  gradient  is  simply  the  ve¬ 
locity,  i.e.  V<^=v. 

Following  Liggett  and  Liu,‘^  the  BEM  formulation  of 
Laplace’s  equation  becomes 


a<t>{ri)  + 


r  dG 

Jr  dn 


■qG 


dr  =  0; 


(1) 


where  0(r)  is  the  value  of  the  potential  at  a  point  r„  F 
denotes  the  boundary  of  the  domain,  G  is  the  Green’s  func¬ 
tion  corresponding  to  Laplace’s  equation,  n  is  the  outward 
normal  of  the  boundary  and  q  =  d<l>ldn  is  the  velocity  normal 
to  the  surface.  The  quantity  a  results  from  singular  contri¬ 
butions  due  to  an  integration  over  the  “base  point”  and  is 
equal  to  tt  for  nodes  along  a  smooth  boundary.  At  a  sharp 
comer,  a  is  equal  to  the  angle  of  the  comer  as  measured 
from  inside  the  domain.’ 
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FIG.  1.  Computational  domain. 


The  free  space  Green’s  function  solution  to  the  axisym- 


metric  Laplacian  is 
^rK(p) 


a  =  (r+r,)-  +  (z-c,)- 


P  “  /  ■  .  7 - n 

where  K{p)  is  the  complete  elliptic  integral  of  the  first  kind. 
The  quantity  SGIdn  can  be  expressed: 


dG  -2 

—  =  —  nMp)  + 


(r“r,r  +  (z-Z,r 


X[(r2-r. -(z-z,)^)n^+2rn,(z-z,)]j,  (5) 

where  E{p)  is  the  complete  elliptic  integral  of  the  second 
kind. 

Since  Eq,  (1)  involves  an  integration  only  around  the 
boundary,  we  need  not  discretize  the  entire  domain.  To  per¬ 
form  this  integration,  it  is  necessary  to  assume  a  behavior  of 
<j>  and  q  over  the  length  of  an  element,  and  we  assume  a 
linear  variation  between  “nodes”  in  our  model.  By  assum¬ 
ing  a  linear  variation  between  nodes,  the  integration  in  Eq. 
(1)  becomes  a  function  of  geometry  alone,  and  the  governing 
equation  can  be  written  in  matrix  form:^ 

[D]{<^}  =  [5]{?}.  (6) 

where  the  elements  of  5  result  from  the  integration  of  G,  the 
elements  of  D  result  from  the  integration  dG/dn  and  the 
values  of  a  have  been  incorporated  into  the  D  matrix.  Pre¬ 
suming  that  or  ^  is  specified  at  each  node  on  the  boundary, 
this  equation  can  be  solved  for  the  remaining  flow  variable  at 
each  node. 

Figure  1  shows  the  computational  domain  for  the  liquid 
jet.  Nodes  are  fixed  in  the  orifice  with  a  prescribed  inflow 
velocity,  so  that  the  BEM  solution  returns  the  velocity  po¬ 


tential.  The  velocity  potential  is  specified  along  the  free  sur¬ 
face  of  the  jet,  and  q  is  returned  by  the  BEM  solver.  The 
surface  is  tracked  by  allowing  the  nodes  to  move  with  their 
local  velocity  vectors.  In  this  case,  flow  kinematics  require 

Dz  d<t>  Dr  d<t> 

Here.  d<i>/dz  and  d<t>fdr  can  be  written  in  terms  of  deriva¬ 
tives  parallel  and  normal  to  the  surface  (d(l>fds  and  q)  and 
the  local  wave  slope  through  a  standard  coordinate  transfor¬ 
mation.  The  derivative  d<t>ISs  is  obtained  via  a  fourth-order 
centered  difference  approximation.  The  unsteady  Bernoulli 
equation  provides  the  boundary  conditions  along  a  free  sur¬ 
face  interface.  Using  the  liquid  density  (p),  orifice  radius 
(a)  and  mean  inflow  velocity  {U)  as  dimensions,  the  non- 
dimensional  form  of  Bernoulli’s  equation  can  be  written: 

d(t>  \  .  K  Bo 

^  =  (8) 

where  We- pU^alcr  is  the  Weber  number.  Bo- pga‘f(T  is 
the  Bond  number,  <j  is  the  liquid  surface  tension  and  k  is  the 
local  surface  curvature.  The  Eulerian-Lagrangian  transfor¬ 
mation  for  nodes  on  the  interface  moving  with  the  liquid 
velocity  is 

DO  dO 

—  =  -  +  V0.V().  (9) 

Using  this  transformation,  the  Bernoulli  equation  in  the  liq¬ 
uid  becomes 

D<f>  I  .  K  Bo 

(10) 

Equations  (7)  and  (10)  are  integrated  in  time  using  a  4th- 
order  Runge-Kutta  scheme  to  solve  for  the  evolution  of  the 
liquid  jet.  To  keep  a  roughly  constant  grid  spacing,  cubic 
splines  are  used  to  regrid  the  surface  of  the  jet  every  time 
step. 

As  the  jet  issues  from  the  orifice,  a  droplet  will  be 
formed  at  the  tip  of  the  jet  and  the  radial  coordinate  of  the 
necking  region  will  tend  to  zero.  To  extend  calculations  be¬ 
yond  droplet  pinching  events,  a  “pinch  criteria”  is  required 
to  remove  droplets  from  the  main  body  of  the  jet  at  times 
when  the  radial  coordinate  in  the  necking  region  is  still  a 
finite  positive  value.  In  the  following  calculations,  we  pre¬ 
sume  that  droplet  pinching  occurs  when  the  minimum  radius 
in  the  necking  region  is  less  than  5%  of  the  orifice  radius. 
When  this  event  occurs,  the  node  at  this  location  is  moved  to 
the  symmetry  axis,  and  the  value  of  the  velocity  potential  is 
estimated  based  on  a  zeroth-order  extrapolation.  Approaches 
similar  to  this  have  been  verified  by  other  researchers.^’ In 
addition,  numerical  experiments  utilizing  differing  values  of 
the  pinch  criteria  have  verified  that  near  the  selected  value, 
droplet  sizes  and  velocities  are  insensitive  to  this  parameter. 

III.  MODEL  VALIDATION 

The  axisymmetric  BEM  solver  and  free  surface  modules 
have  been  validated  with  results  for  the  oscillation  of  a  liquid 
droplet  in  the  absence  of  gravity.^  In  the  linear  regime,  os- 
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FIG.  2,  Jet  profile  comparison  near  the  pinch  location,  (a)  experiment  of 
.Moses,  lb)  BE.M  Result.  lVe=l7.6.  Bo=0.0l09.  *  =  0.447. 


cillation  frequencies  are  within  one  tenth  of  one  percent  of 
the  analytical  frequency  given  by  Lamb.'"*  For  moderate  os¬ 
cillations.  the  BEM  solver  calculates  frequency  shifts  rela¬ 
tive  to  the  linear  frequency  within  0.5%  of  the  analytic  result 
presented  by  Tsamopoulos  and  Brown. 

Moses'  investigated  the  breakup  of  a  liquid  jet  with  an 
experimental  setup  using  water  in  air.  He  incorporated  digital 
image  processing  techniques  to  provide  accurate  data  on 
droplet  size  and  surface  wave  shapes.  His  setup  consisted  of 
a  tank  draining  through  an  orifice  with  radius  a  =  0.0283  cm 
and  length  to  orifice  radius  ratio  Ua  =  10.22.  A  speaker  ex¬ 
ternally  fitted  to  the  tank  with  an  input  sine  wave  from  a 
signal  generator  provided  acoustic  perturbations  to  the  jet. 
The  jet  was  illuminated  with  a  spark  source,  and  imaged  by 
a  CCD  camera  connected  to  a  video  tape  recorder.  Edge 
detection  was  by  the  convolution  of  a  video  image  of  the  jet 
with  a  mask.  Data  were  collected  at  fixed  Weber  numbers  of 
17.6  and  76.6,  and  with  wave  numbers  in  the  range 
0.3 0.85.  Figure  2(a)  shows  a  typical  jet  image  near  the 
pinch  location  with  satellite  droplet  formation. 

A  BEM  calculation  was  performed  in  order  to  match  the 
experimental  results  of  Moses.  The  numerical  simulation 
used  a  grid  spacing  of  20%  of  the  orifice  radius,  a  time  step 
Af  =  0.0025  and  separated  droplets  when  the  radial  coordi¬ 
nate  at  the  pinch  point  was  within  5%  of  the  orifice  radius. 
Figure  2(b)  shows  a  three-dimensional  rendering  of  the  jet 
profile  at  time  f  =  92.4  time  units.  Here,  the  BEM  result  was 
aligned  axially  with  the  experimental  image  such  that  the 
breakup  point  occurs  at  the  same  axial  location;  the  inviscid 
simulation  tends  to  underpredict  jet  breakup  length  due  to  the 
neglect  of  viscous  effects.  This  figure  shows  the  excellent 
qualitative  agreement  between  the  experimental  and  numeri¬ 
cal  results.  Note  how  the  tip  of  the  jet  shows  a  series  of  three 
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FIG.  3.  Typical  jet  evolution.  We  =  l00.  *  =  0.70.<7' =2%. 


or  four  ripples  formed  by  the  recoil  of  the  surface  after  the 
pinching  event.  In  our  simulation,  these  disturbances  lead  to 
atomization  of  the  jet  tip  into  3-4  sub-satellite  droplets.  Inte¬ 
gration  of  the  motion  of  these  sub-satellite  droplets  show 
their  tendency  toward  recombination.  That  is,  the  satellite 
droplet  shown  in  Figure  2(b)  is  actually  the  superposition  of 
two  sub-satellite  droplets  which  were  integrated  in  time 
separately,  and  therefore  the  size  of  droplet  shown  is  less 
than  the  actual  size  due  to  this  overlap.  The  model  is  cur¬ 
rently  incapable  of  recombining  droplets  after  pinch. 


To  simulate  experimental  results  of  liquid  jet  breakup, 
the  computational  domain  shown  in  Figure  I  was  used  with 
an  inflow  velocity  along  the  entire  orifice  exit  of  the  form 

(7  =  (/[l+,?'sin(kr)],  (H) 

where  q'  is  the  amplitude  of  the  oscillatory  disturbance  and 
k  is  the  wave  number  of  the  perturbation.  Calculations  were 
run  for  various  velocity  perturbations,  wave  numbers,  and 
Weber  numbers.  Gravity  is  neglected,  so  that  Bo  =  Q.  All 
calculations  begin  with  a  long  column  of  fluid,  usually  about 
15  to  20  orifice  radii,  outside  the  orifice  in  order  to  ensure 
that  perturbations  from  the  assumed  spherical  end  cap  do  not 
influence  the  development  of  waves  due  to  the  imposed  per¬ 
turbation.  The  calculations  presented  here  use  a  grid  spacing 
of  20%  of  orifice  radius,  a  time  step  of  At =0.0025  time 
units  (unless  noted  otherwise),  and  separate  droplets  from  the 
calculation  when  the  radial  coordinate  at  the  pinch  location  is 
within  5%  of  the  orifice  radius. 

A  standard  case  is  chosen  as  We  =100,  k  =  0.7  and 
q' =  2%.  Linear  stability  analysis  shows  that  this  wave  num¬ 
ber  is  the  most  unstable  for  this  Weber  number.’’  Figure  3 
shows  jet  profiles  for  five  successive  droplet  shedding  events 
for  the  standard  case,  beginning  at  f =90.0  time  units.  By  this 
time,  the  calculated  drop  sizes  and  jet  breakup  length  pre¬ 
dicted  by  the  BEM  solver  have  a  repeating  pattern  with  very 
little  variation.  Periodic  bulges  appear  in  the  jets  due  to  the 
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TABLE  1.  Jet  breakup  simulation  results. 


We 

k 

50 

0.7 

2% 

1.845 

0.713 

100 

0.5 

1.946 

1.250 

100 

0.7 

2% 

1.818 

0.867 

100 

0.9 

2% 

1.719 

0.443 

100 

l.l 

2% 

1.616 

0.183 

100 

0.7 

1.793 

0.986 

100 

0.7 

4^ 

1.771 

1.053 

100 

0.7 

5^ 

1.762 

1.091 

100 

0.7 

6^f 

1.750 

1.136 

100 

0.7 

l^c 

1.749 

1.158 

unsteady  inflow,  as  noted  by  Reba  and  Brosilow.‘®  The  non¬ 
linear  effects  on  the  bulges  lead  to  the  formation  of  “main" 
and  “satellite"  droplets,  a  behavior  well  known  for  low- 
speed  jets. 

A  moderate  29c  disturbance  magnitude  was  chosen  over 
typical  acoustically-generated  perturbation  magnitudes  in  or¬ 
der  to  keep  the  jet  breakup  length  to  a  reasonable  size,  and 
therefore  make  calculation  times  acceptable.  After  the  jet 
sheds  the  fluid  from  the  initial  assumed  geometry,  there  were 
between  255  and  290  nodes  on  the  jet  surface  during  this 
calculation.  The  first  77  time  units  of  the  calculation  required 
about  500,000  CPU  seconds  on  an  IBM  RISC/6000  com¬ 
puter.  Table  I  presents  a  summary  of  the  droplet  radii  for  the 
calculations  in  this  report. 

A.  Effect  of  oscillation  magnitude 

The  effect  of  the  oscillation  magnitude  on  the  jet 
breakup  is  investigated  by  six  simulations  of  q'  =  2% 
through  7%.  Figure  4  shows  jet  profiles  for  three  of  these 
disturbance  magnitudes.  Each  pair  of  profiles  shows  the  jet  at 


(a) 


(c) 


FIG.  4.  Effect  of  inflow  velocity  disturbance  magnitude,  W^=100, 
je  =  0.70,  (a)  (7'=2%,  (b)  ^'  =  4%,  (c)  q'  =  6%. 


RG.  5.  Droplet  radii  vs.  inflow  velocity  disturbance  magnitude  for  We 
-  100.  k  =  0.1  (o  finite-length  calculation: _ infinite  jet  calculation). 


times  just  before  a  main  droplet  and  a  satellite  droplet  are 
shed  from  the  calculation.  As  this  figure  shows,  increasing 
the  size  of  the  perturbation  decreases  the  breakup  length  and 
changes  the  shape  of  the  droplets.  At  larger  disturbance  mag¬ 
nitudes,  the  main  droplets  take  on  a  “squashed"  shape  as  a 
result  of  the  high  velocity  fluid  (when  -Trl2+2mr 
<  kt<rr/2+2mr)  encountering  lower  velocity  fluid  which 
has  already  exited  the  nozzle.  This  phenomena,  known  as  the 
“Klystron"  effect,  has  been  documented  qualitatively  by  nu¬ 
merous  researchers  in  the  case  where  forcing  perturba¬ 
tions  are  very  large  amplitude. 

In  the  other  half  of  the  imposed  oscillation  (when 
7r/2  +  2n7r<^r<3Tr/2+2n7r),  the  jet  is  “stretched"  since 
the  fluid  immediately  outside  the  orifice  is  at  a  higher  veloc¬ 
ity  than  the  inflow.  Therefore,  the  overall  effect  of  increased 
disturbance  amplitude  is  to  narrow  the  peaks  and  broaden  the 
troughs  in  the  surface  waves  on  the  jet.  Breakup  times  are 
reduced  (as  compared  to  the  infinitesimal  case)  primarily  due 
to  the  fact  that  the  troughs  begin  at  a  lower  height. 

Calculated  droplet  sizes  are  compared  with  an  infinitesi¬ 
mal  disturbance^  values  in  Figure  5.  The  satellite  size  is 
shown  to  increase  with  disturbance  magnitude  at  the  wave 
number  ratio  selected  for  these  simulations.  Necking  and 
pinching  of  the  jet  under  finite-amplitude  disturbance  occur 
at  locations  much  closer  to  the  peak  of  the  wave.  The  broad, 
flat  troughs  created  by  these  perturbations  explain  the  in¬ 
crease  in  satellite  droplet  sizes  with  increasing  disturbance 
magnitude. 

Extrapolation  to  higher  values  would  presumably 
lead  to  a  monodisperse  case  in  which  both  drops  are  the 
same  size.  Through  the  use  of  piezoelectric  drivers,  droplet 
trains  of  this  type  have  been  created  experimentally  by 
Dressier  and  coworkers.^  Direct  comparison  with  these  ex¬ 
periments  is  not  possible  because  pressures  are  not  typically 
measured  in  the  plenum  upstream  of  the  orifice  in  Dressier' s 
device. 
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FIG.  6.  Effect  of  disturbance  wave  number.  We  =100,  q'  =  2%,  (a) 
4  =  0.5.  (b)  4  =  0.7,  (c)  4  =  0.9,  (d)  4=  1.1.  (e)  4=  1.3. 


RG.  7.  Droplet  radii  vs.  disturbance  wave  number. 


B.  Effect  of  oscillation  wave  number 

Since  the  jet  radius  and  mean  orifice  exit  velocity  are 
chosen  as  characteristic  dimensions,  non-dimensional  wave 
numbers  are  equivalent  to  disturbance  frequency.  Figure  6 
shows  the  effect  of  wave  number  on  the  character  of  the  jet 
for  We  =  100  and  q*  =  2%,  Each  of  the  first  four  pairs  of 
profiles  show  the  jet  at  times  just  before  a  main  droplet  and 
just  before  a  satellite  droplet  are  shed  from  the  calculation. 
Increasing  the  wave  number  from  ^  =  0.5  to  1.1  tends  to 
decrease  the  breakup  length,  even  beyond  the  most  unstable 
wave  number,  ^,„ax  =  0.7;  a  trend  not  predicted  by  linear 
theory.  In  addition,  the  size  of  satellite  drops  tends  to  de¬ 
crease  with  increasing  frequency  in  this  range.  At  a  wave 
number  of  /:=  l.l,  satellite  drops  have  nearly  vanished,  in¬ 
dicating  that  a  wave  number  near  this  value  can  produce 
monodisperse  atomization.  In  fact,  viscous  effects  would 
probably  inhibit  the  formation  of  the  very  small  satellite 
droplet  in  Figure  6(d).  At  /:=  1.3,  the  breakup  length  of  the 
liquid  jet  increases  drastically  from  the  previous  calculation, 
and  the  model  no  longer  predicts  the  formation  of  satellite 
droplets.  Figure  6(e)  shows  the  jet  profile  at  time  r=  107.0 
time  units;  the  jet  length  length  is  still  increasing  at  this  point 
in  the  calculation.  Note  the  damping  of  the  waves  along  the 
jet  surface. 

Quantitative  predictions  of  droplet  sizes  are  shown  in 
Figure  7.  As  in  Figure  5,  the  solid  lines  are  results  for  infini¬ 
tesimal  disturbances.  Here,  experimental  results  are  a  com¬ 
posite  of  data  from  Lafrance^  and  Rutland  and  Jameson.^ 
The  finite-amplitude  results  are  for  <7'  =  2%,  and  they  show 
a  tendency  to  drive  conditions  toward  a  monodisperse  result 
by  increasing  satellite  drop  size  and  decreasing  main  drop 
size  over  the  entire  range  of  wave  numbers  investigated. 
Once  again,  the  peak-sharpening  and  trough-broadening 
which  occur  with  increased  disturbance  magnitude  provide 
an  explanation  for  the  observed  results.  At  higher  distur¬ 
bance  amplitudes,  local  curvature  in  the  transition  region  be¬ 
tween  peaks  and  troughs  can  be  driven  to  large  enough  val¬ 


ues  to  promote  nonlinear  instability  at  wavelengths  predicted 
to  be  stable  using  linear  analysis. 

C.  Effect  of  Weber  number 

For  fixed  orifice  radius  and  liquid  density,  decreasing  the 
Weber  number  corresponds  to  either  increasing  the  surface 
tension  or  decreasing  the  inflow  velocity.  As  the  surface  ten¬ 
sion  is  increased,  the  wave  formed  on  the  surface  of  the  jet 
from  the  unsteady  inflow  should  tend  to  grow  more  quickly, 
leading  to  shorter  breakup  lengths.  This  effect  is  shown  in 
Figure  8  for  We  =  50  and  100,  using  k  =  0.7  and  ^'  =  2%. 
Each  pair  of  profiles  shows  the  jet  just  before  a  main  droplet 
and  satellite  droplet  are  shed  from  the  calculation.  Increasing 
the  Weber  number  also  tends  to  decrease  the  size  of  the 
satellite  droplets  for  the  chosen  conditions.  Note  that  the  tip 
of  the  satellite  droplet  for  the  lower  Weber  number  case 
takes  on  a  more  spherical  shape  before  being  shed  from  the 
calculation.  In  addition,  at  lower  We  values,  the  pinching 
event  tends  to  form  a  series  of  waves  on  the  parent  surface 
due  to  the  enhanced  influence  of  surface  tension.  This  effect 


(a) 


(b) 


RG.  8.  Effect  of  Weber  number  on  jet  profile.  4  =  0.7,^' =  2%,  (a)  We 
=  100,  (b)  We  =50. 


1578  Phys.  Fluids,  Vol.  8,  No.  6,  June  1996 


J.  H.  Hilbing  and  S.  0.  Heister 


28 


1 


I 


Satellite  Droplet 

time  =  69.6  Rebound  OistafX» 


no.  9.  Satellite  droplet  evolution  just  prior  to  shedding. 


can  also  be  observed  (both  in  calculations^  and  experimental 
results*)  in  the  low  We  case  shown  in  Fig.  2. 

D.  Droplet  velocities 

Infinite  jet  simulations  lead  to  droplets  with  no  velocity 
relative  to  each  other  due  to  the  assumption  of  periodicity 
over  a  given  wavelength  of  the  jet.^  However,  for  finite- 
length  liquid  jets,  droplets  can  be  shed  sequentially,  so  that 
the  droplet  velocities  are  not  necessarily  equal  to  each  other 
or  to  the  orifice  velocity.  The  calculations  presented  here 
show  that  the  satellite  droplets  are  shed  from  the  jet  with 
velocities  less  than  that  of  the  main  droplets.  Figure  9  shows 
the  jet  profiles  at  three  times  near  a  droplet  pinch.  Each  pro¬ 
file  in  this  figure  has  been  shifted  axially  so  that  the  left-hand 
side  of  the  jet  remains  in  the  same  relative  location.  Just  after 
a  main  droplet  is  shed  from  the  jet,  the  surface  of  the  satellite 
droplet  is  roughly  conical  at  the  tip.  Surface  tension  acts  to 
pull  this  liquid  into  a  spherical  shape,  and  tends  to  slow  the 
velocity  of  the  satellite  droplet.  Note  in  Figure  9  at  r  =  69.6 
time  units  the  distance  that  the  surface  has  rebounded  to¬ 
wards  the  orifice  before  the  satellite  droplet  is  separated  from 
the  calculation.  This  phenomenon  is  also  present  for  the 
main  droplets,  so  that  a  monodisperse  droplet  train  should 
have  droplets  of  constant  axial  velocity  slightly  less  than  the 
mean  orifice  velocity. 

The  droplet  velocity  difference  is  defined  as  the  percent 
difference  in  satellite  and  main  droplet  velocities  relative  to 
the  main  droplet  velocity.  This  difference  is  calculated  by 
integrating  the  motion  of  the  droplet  after  being  shed  from 
the  jet.  The  position  of  the  droplet’s  calculated  center  of 
mass  is  differenced  to  determine  the  droplet  velocity.  Figure 
10  shows  a  typical  plot  of  the  position  and  velocity  of  the 
center  of  mass  for  a  main  droplet  and  a  satellite  droplet. 
Table  II  presents  velocity  differences  for  selected  calcula- 
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no.  10.  Droplet  center  of  mass  position  and  velocity  iWe  =  100. 
ifc=0.7.^'  =  2^). 


tions.  As  this  table  shows,  the  velocity  difference  increases 
when  the  Weber  number  decreases,  the  wave  number  in¬ 
creases  or  the  disturbance  magnitude  decreases.  This  in¬ 
crease  in  velocity  difference  can  be  attributed  to  increasing 
surface  tension  or  decreasing  satellite  droplet  size, 

Pimbley  and  Lee**  investigated  satellite  droplet  forma¬ 
tion  experimentally,  noting  the  velocity  difference  between 
the  main  and  satellite  droplets.  They  found  appropriate  exci¬ 
tation  conditions  to  separate  satellite  droplets  from  either 
side  of  a  main  drop,  as  well  as  forming  both  droplets  simul¬ 
taneously,  the  so-called  “infinite”  satellite  condition.  Satel¬ 
lite  droplets  not  at  the  infinite  condition  were  found  to  merge 
with  main  droplets  a  few  wavelengths  downstream.  At  the 
lowest  excitation  amplitudes,  rear-merging  satellites  were 
observed  in  the  experiment.  These  satellite  droplets  merged 
backwards  with  next  main  droplet.  As  the  amplitude  was 
increased,  the  velocity  difference  decreased  until  the  infinite 
condition  was  reached.  Increasing  the  amplitude  farther  re¬ 
sulted  in  forward-merging  satellites,  which  merged  with  the 
parent  main  droplet.  The  BEM  calculations  presented  here 
agree  with  the  low-amplitude  experimental  results  which 
produced  backward-merging  satellite  droplets. 


E.  Amplitude  modulated  disturbances 

Orme  and  Muntz  studied  the  breakup  of  a  liquid  jet 
using  amplitude-modulated  sinusoidal  disturbances.  They 


TABLE  II.  Droplet  velocity  differences. 


We 

k 

Velocity  difference 

50 

0.7 

2% 

-12.0±1.4% 

100 

0.5 

2% 

-l.4i0.9% 

iOO 

0.7 

2% 

-3.4±0.9% 

100 

0.9 

2% 

-3.9  ±0.9% 

100 

1.1 

2% 

-9.8±4.7% 

100 

0.7 

4% 

-l.9±0.6% 
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no.  II.  Jet  breakup  with  an  amplitude-modulated  inflow  velocity,  (a)  in¬ 
flow  disturbance  waveform,  (b)  jet  profiles  at  four  successive  pinch  events, 
(c)  relative  droplet  velocities. 


found  that  the  droplets  coalesced  downstream  due  to  varia¬ 
tions  in  their  velocities.  For  this  type  of  disturbance,  the 
inflow  velocity  can  be  written  in  the  form 

<?  =  i/|n-  ^'  +  9;sin|^|  sin(/cr)  ’,  (12) 

where  is  the  magnitude  of  the  modulated  disturbance  and 
is  the  modulation  frequency  ratio. 

A  BEM  calculation  was  performed  in  order  to  simulate 
these  experiments.  We  chose  disturbance  magnitudes  of 
<7'  =  0.015  and  =  0.0074,  which  correspond  to  a  5%  fluc¬ 
tuation  in  the  reservoir  pressure,  a  wave  number  of  it  =  0.7 
and  a  frequency  ratio  of  A^=4.  Figure  1 1(a)  shows  the  inflow 
disturbance  waveform  corresponding  the  these  values.  A 
Weber  number  of  W^  =  200  is  chosen  for  this  simulation  in 
order  to  limit  the  CPU  time  required;  gravity  is  neglected. 
The  jet  begins  with  length  1 5a  outside  the  orifice,  and  is 
marched  with  a  time  step  of  Ar=0.01  through  a  number  of 
droplet  shedding  events. 

Orme  and  Muntz  predict  the  shedding  of  a  fluid  packet 
containing  JV  pairs  of  main  and  satellite  droplets,  the  subse¬ 
quent  breakup  of  this  packet  into  2N  droplets,  followed  by 
the  recombination  of  these  droplets  into  one  large  droplet. 
Figure  11(b)  shows  the  jet  profiles  at  four  times  just  prior  to 
four  successive  pinch  events.  These  four  pinch  events  in¬ 
clude  the  eight  droplets  for  one  complete  cycle  of  the  modu¬ 
lation  frequency.  Instead  of  a  single  packet  with  eight  drop¬ 
lets,  we  find  that  three  packets  of  one  droplet,  plus  one 
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packet,  at  r  =  1 3 1 .0  time  units,  of  two  main  and  three  satellite 
droplets.  When  this  multiple  droplet  packet  is  integrated  after 
pinch,  the  two  droplets  on  the  orifice-side  of  the  packet  pinch 
off.  The  velocity  of  each  droplet's  center  of  mass  is  calcu¬ 
lated,  and  the  arrows  in  Figure  1 1(c)  show  their  direction  of 
motion  relative  to  the  fluid  packet  consisting  of  one  main  and 
two  satellite  droplets.  Our  results  of  Fig.  1 1(c)  compare  well 
with  Fig.  2  from  Orme  and  Muntz,  predicting  the.  formation 
of  a  single  droplet. 


V.  CONCLUSIONS 

A  numerical  model  has  been  developed  to  investigate  the 
influence  of  unsteady  inflow  conditions  on  the  nonlinear  evo¬ 
lution  and  droplet  formation  processes  within  a  low  speed, 
finite-length  liquid  jet.  The  size  of  both  main  and  satellite 
droplets  is  predicted  under  the  assumption  of  a  sinusoidal 
perturbation  to  the  orifice  exit  velocity.  Results  indicate  that 
modulation  of  either  the  amplitude  or  the  frequency  (wave 
number)  of  the  perturbation  can  affect  droplet  sizes  so  as  to 
create  a  monodisperse  droplet  train.  In  addition,  velocity  dis¬ 
persions  (between  main  and  satellite  droplets)  are  predicted 
with  the  model;  therefore,  the  tendency  of  droplets  to  recom¬ 
bine  after  pinch  off  can  also  be  addressed. 
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8  Appendix  C  -  Acoustic  Interactions  with  Liquid  Jet 
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Nomenclature 

a  =:  orifice  radius 

Im  =  modified  Bessel  fct.  of  1st  kind,  order  m 
Km  =  modified  Bessel  fct.  of  2nd  kind,  order  m 
k  =  wave  number 

n  =  mode  number  for  column  excitation  (Eq.  5) 

P  =  pressure 

q  =  velocity  normal  to  local  surface 
r  =  radial  coordinate 
t  =  time 
U  =  velocity 

We  =  Weber  number,  We  —  pgU'^afcr  =  pU^afa  in  Figs.  8-12 

X  =  horizontal  coordinate 

y  =  vertical  coordinate 

z  =  axial  coordinate 

K  =  surface  curvature 

p  =  density 

e  —  density  ratio,  pg/ p 

cr  =  surface  tension 

Lj  =  oscillation  frequency 

(j)  zz  velocity  potential 

Subscripts 
0^  gas  phase 

Superscripts 

0  =:  perturbation  quantity 
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Abstract 


The  effect  of  unsteady  chamber  characteristics  on  the  atomization  process  is  investigated  using  non¬ 
linear  free-surface  models  based  on  Boundary  Element  Methods  (BEMs).  Two  separate  scenarios 
are  considered.  In  the  first  case,  a  liquid  column  is  subjected  to  transverse  oscillations  from  the  gas 
phase.  The  second  calculation  addresses  ‘'dynamic  orifice  massflow”  due  to  the  presence  of  longitu¬ 
dinal  acoustic  perturbations.  Results  indicate  that  the  atomization  process  can  be  strongly  affected 
by  these  perturbations.  As  expected,  the  most  prominent  coupling  occurs  at  driving  frequencies 
which  are  at  (or  near)  the  natural  frequency  of  oscillation  of  the  liquid. 


Introduction 

In  a  liquid  rocket  engine  (LRE)  combustion  chamber,  the  atomization  process  serves  as  a  precursor 
to  complex  vaporization,  mixing,  and  reaction  processes.  Not  only  is  the  atomization  process  im¬ 
portant  to  characterizing  the  steady-state  performance  of  a  LRE,  but  it  also  can  play  an  important 
role  in  unsteady  processes  within  the  combustion  chamber.  In  fact,  numerous  authors^“^  have 
implicated  the  atomization  process  as  a  mechanism  to  (at  least  partially)  explain  high-frequency 
LRE  combustion  instabilities. 

Even  under  steady  conditions,  our  knowledge  of  the  detailed  processes  leading  to  atomization  of 
a  liquid  jet  is  quite  modest.  For  this  reason,  there  have  been  relatively  few  efforts  aimed  at  improving 
our  understanding  of  injection  processes  under  dynamic  conditions.  Most  previous  efforts  have  been 
experimental  in  nature  and  have  been  motivated  by  LRE  combustion  stability  concerns. 

One  group  of  experiments  have  focused  on  longitudinal  acoustic  fields^*"^  investigating  both 
low  injection  velocities'’^  and  high  injection  velocities®’^  characteristic  of  actual  LRE  injectors. 
These  authors  concluded  that  the  presence  of  a  time-dependent  pressure  field  at  the  injection  point 
leads  to  a  dynamic  orifice  massflow,  even  in  the  event  that  the  injector  manifold  pressure  remains 
constant.  At  low  injection  velocities,  this  effect  leads  to  periodic  “bulges”  in  the  jet  diameter;  this 
character  also  persists  into  the  higher  jet  velocity  regime®.  In  addition,  droplet  coalescence  was 
effected  by  the  time-dependent  injection  flowrate.  Similar  effects  have  been  observed^®’^^  when  the 
chamber  pressure  is  held  fixed,  but  the  flowrate  is  varied  using  piezoelectric  drivers  upstream  of 
the  orifice  inlet.  At  high  injection  velocities,  Ingebo^  showed  a  dramatic  decrease  in  mean  droplet 
size  with  the  amplitude  of  the  disturbance.  This  compelling  evidence  indicates  the  importance  of 
the  acoustic  interaction  with  the  atomization  process. 

A  second  group  of  researchers®’^^"^"^  have  focused  on  interactions  of  a  liquid  jet  with  a  transverse 
acoustic  perturbation.  Here,  experiments  have  revealed  a  broadening  of  the  jet  cross-section^^,  and 
a  physical  deflection  of  the  jet  due  to  the  “crosswind”  from  the  acoustic  field^®.  Here,  the  variable 
massflow  effect  can  also  be  present  due  to  the  fact  that  the  injector  manifold  reacts  to  the  average 
pressure  at  the  injector  face  and  distributes  fluid  spatially  in  accordance  with  the  local  injector  face 
pressure^^. 

Analytic  models  aimed  at  quantifying  the  phenomena  described  above  have  made  heavy  use  of 
experimental  results.  Early  efforts®""^  made  use  of  one-dimensional  mass  and  momentum  balances 
in  which  liquid  and  gas  phases  were  decoupled.  Other  researchers  have  simply  sought  to  correlate 
their  experimental  results  with  injection  and  acoustic  wave  conditions.  More  recently,  the  jet 
broadening  phenomena  present  under  some  transverse  acoustic  interactions  was  quantified  using  a 
static  equilibrium  analysis^"^. 

While  these  efforts  have  advanced  our  understanding  of  these  processes,  the  effect  of  coupling  of 
gas  and  liquid  flowfields  in  a  dynamic  environment  has  yet  to  be  considered.  Recent  developments 
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of  numerical  models  based  on  the  use  of  Boundary  Element  Methods  (BEMs)^'^’^^  permit  us  to 
analyze  the  acoustically-driven  interactions  described  above.  The  use  of  the  BEM  approach  provides 
high  resolution  of  the  interface  (under  very  large  distortions),  as  well  as  the  capability  to  extend 
simulations  beyond  atomization  events.  This  paper  will  present  simulations  of  both  longitudinal 
and  transverse  acoustic  wave  interactions  with  an  initially-cylindrical  column  of  fluid.  Transverse 
simulations  will  include  coupling  of  gas  and  liquid  velocity  fields  and  longitudinal  simulations  will 
include  the  effect  of  dynamic  orifice  massflow  in  time-accurate  calculations  valid  for  nonlinear 
deformations  of  the  liquid. 


Classification  of  Acoustically- Driven  Instabilities 


Before  performing  the  numerical  simulations  described  above,  it  is  prudent  to  determine  the  applica¬ 
ble  flow  regimes  in  LREs,  as  well  as  those  studied  by  other  researchers.  One  important  consideration 
involves  the  “frequency  response”  of  the  liquid  to  the  imposed  oscillation.  Consider  the  cylindrical 
column  of  liquid  subject  to  either  longitudinal  or  transverse  acoustic  perturbations  as  shown  in  Fig. 
1.  Variables  shown  in  this  figure  are  defined  in  the  Nomenclature. 

For  the  transverse  oscillation  in  Fig.  1,  fundamental  frequencies  are  obtained  from  a  linearized, 
inviscid,  2-D  analysis  similar  to  the  classic  axisymmetric  analysis  used  by  Lamb^^  in  the  case  of  a 
droplet.  In  the  2-D  case,  it  is  straightforward  to  show  that  velocity  potentials  satisfying  Laplace’s 
equation  in  the  liquid  and  gas  phases  can  be  represented: 


<j>  =  — (r/a)"^  cos(n^)  cos(u;^)  (1) 

CiU) 

4>g  =  ——{a/r)^cos{nO)cos{u}t)  (2) 

assuming  a  surface  shape  r  =  a  +  dcos(n^)  sin(tj<).  Here,  n  is  the  order  of  the  oscillation  and  we 
will  assume  that  d  is  a  perturbation  (d  <<  a).  Since  the  column  is  initially  at  rest,  the  linearized 
result  from  Bernoulli’s  equation  can  be  expressed: 


^  a  ^  dt  dt 


=  K<J 


where  k  is  the  linearized  curvature  of  the  distorted  surface: 


(3) 


K  =  i[l  ~  ~(1  —  n^)  cos(n9)  sin(c<;^)] 
a  a 

Combining  Eqs.  1-4,  we  obtain  the  fundamental  frequencies  of  oscillation  of  the  column: 


(4) 


n  n(n‘^  —  l)cr 


{P  +  PgW 


(5) 


As  in  the  case  of  a  droplet^^,  we  see  that  the  lowest  order  mode  is  for  n  =  2.  In  this  Ccise,  the 
fundamental  frequency  is  given  by  =  Q<7/{{p  +  Pg)a^)  which  gives  a  result  13%  lower  than  that 
of  a  droplet  in  a  low  density  gas. 

For  the  case  of  a  longitudinal  excitation,  the  fundamental  frequency  of  the  liquid  column  can  be 
determined  by  analysis  of  a  dispersion  relation^®  which  describes  the  variation  of  u  with  changes 
in  wave  number,  k: 


,  lojk)  Kojk) 

iu{k) 


)  +  ‘2iku}  —  k 


J<o{k)  k{l-k^) 
Ki{k)  ~  We 


=  0 


(6) 
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Table  1:  Frequency  Characterization,  Acoustic  Interactions 


Researcher(s) 

AP/P 

cxlO^ 

ire 

Frequen 

^9 

cies  in  Hz 

Reba  h  Brosilow^  (L) 

0.03-0.49 

1.2-4.2 

0.38-21 

250-500 

4-2700 

Ingebo^  (L) 

0.11-0.22 

0.5-4. 4 

47-1700 

1190 

600-373,000 

Miesse"*  (T) 

— 

1.2 

— 

1100-8800 

100-400 

Hoover,  et.  (T) 

0.0046-0.055 

1.2-14 

0.03-80 

500-2500 

470-920 

Oefelein  k.  Yang^^  (T) 

0.65-4.0 

17.8 

— 

454-538 

80  -  270* 

Note:  (T),  (L)  denote  tangential  and  longitudinal  modes,  respectively, 
*  -  Fuel  and  oxid.  orifices  for  5U  and  Double-Row  Cluster  injectors. 


Here,  we  assume  that  variables  have  been  nondimensionalized  using  the  jet  radius,  jet  velocity, 
and  liquid  density  as  dimensions.  Under  this  assumption,  the  density  ratio,  e  =  Pg/p,  and  the 
Weber  number,  We  =  pgU^ajcr,  are  the  relevant  dimensionless  groups.  Equation  6  is  developed 
assuming  a  jet  deformation  of  the  form: 


r  =  1  +  (7) 

<2 

such  that  the  real  part  of  cj  corresponds  to  the  “growth  rate”  of  disturbances  with  wave  number  k 
and  the  imaginary  part  of  uj  represents  the  characteristic  frequency  of  disturbances  traveling  along 
the  surface.  Following  Sterling  and  Sleicher^®,  we  presume  that  after  some  finite  time,  only  the 
disturbance  with  maximum  growth  rate  will  be  present.  By  iterating  on  various  k  values,  we  can 
find  the  wave  number  of  the  disturbance  satisfying  this  criteria  and  determine  the  fundamental 
frequency  response  (imaginary  part  of  u;)  corresponding  to  this  condition. 

Using  these  notions,  we  can  characterize  the  frequency  and  flow  regimes  investigated  by  various 
researchers,  as  well  as  those  present  in  engines  which  have  experienced  combustion  instabilities. 
Data  resulting  from  this  procedure  are  summarized  in  Table  1.  This  table  summarizes  experimental 
results  from  cold  and  combusting^’^^  flows  with  imposed  oscillations.  Ingebo^  ran  a  series 

of  combusting  experiments  under  oscillating  chamber  pressure  conditions  (created  using  a  siren  at 
the  chamber  exit),  while  Oefelein  and  Yang^^  report  on  the  extensive  data  available  from  stability 
testing  of  the  F-1  engine.  This  engine  exhibited  tangential-mode  instabilities  in  many  of  the  early 
injector  designs.  In  the  case  of  transverse  oscillations,  Weber  numbers  in  Table  1  are  based  on 
estimated  peak  crossflow  velocities  and  liquid  frequencies  are  calculated  using  the  primary  (n  =  2) 
mode  in  Eq.  5.  For  longitudinal  oscillations,  Weber  numbers  are  based  on  jet  velocity,  and  a;  values 
are  calculated  from  Eq.  6. 

Results  indicate  that  frequencies  observed  in  firings  (or  utilized  in  tests)  are  within  (or  near)  the 
natural  frequency  range  of  the  liquid  columns.  In  particular,  the  observed  instability  frequencies  in 
the  F-1  engine  are  very  near  the  natural  frequency  of  oscillation  for  fuel  jets  used  in  early  (unstable) 
injector  designs  for  this  engine.  Since  u  a  the  natural  frequency  of  the  column  doubles  with 

a  37%  reduction  in  jet  radius  (presumed  to  occur  with  the  shedding  of  drops  from  the  periphery). 
Therefore,  it  will  be  fruitful  to  investigate  jet  response  near  the  condition  u  ^  ujg  to  investigate 
coupling  between  the  atomization  process  and  the  acoustic  modes  from  the  chamber. 
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Modeling  Approach 


Liquid  jet  atomization  problems  represent  significant  challenges  to  modelers  due  to  the  nonlinear 
nature  of  the  free  surface,  and  viscous  and  turbulence  effects.  Current  models  are  unable  to  address 
all  these  physical  processes.  The  present  analysis  is  intended  to  shed  some  light  on  the  gross  response 
of  liquid  jets  to  both  longitudinal  and  transverse  acoustic  oscillations.  In  particular,  we  will  presume 
that  both  liquid  and  gas  phases  are  incompressible  and  inviscid  with  negligible  turbulence  in  the 
flow.  The  lack  of  viscosity  in  the  gas  implies  that  pressure  distributions  in  separated  regions  will 
not  be  resolved.  In  the  liquid,  viscosity  does  not  alter  the  basic  surface  shape  unless  droplet  sizes 
are  of  the  same  order  as  the  boundary  layer  thickness. 

Many  authors  have  proposed  that  the  gas  field  can  be  treated  as  locally  incompressible  in 
analyzing  the  flow  of  the  gas  around  the  liquid  jet.  This  assumption  is  valid  when  the  wavelength 
of  the  acoustic  disturbance  is  much  larger  than  that  which  is  applicable  to  the  liquid.  In  our  case, 
wavelengths  in  the  gas  phase  are  proportional  to  combustion  chamber  dimensions,  whereas,  the 
wavelength  in  the  liquid  jet  is  proportional  to  the  orifice  size.  Since  the  orifice  is  typically  2-3 
orders  of  magnitude  smaller  than  the  applicable  chamber  dimension,  the  incompressible  gas  phase 
assumption  is  prudent. 

In  this  section,  we  will  choose  the  jet  radius  (a),  liquid  density  (p),  and  a  characteristic  velocity 
(C/)  as  dimensions.  Under  this  nondimensionalization,  the  Weber  number,  We  =  pgU^a/cr  and  the 
gas/liquid  density  ratio,  e  =  pg/p  are  the  two  dimensionless  variables  characterizing  these  flows.  For 
an  inviscid  gas  or  liquid  domain,  velocity  potentials  <i>g,(l>  exist  and  must  satisfy  Laplace’s  equation: 

^^(i>  =  V^g  =  0  (8) 


The  unsteady  Bernoulli  equation  provides  conditions  relating  velocity  potentials  at  the  gas/liquid 
interface: 


=  0 


(9) 


where  Pg  is  the  dimensionless  gas  pressure  at  the  interface,  and  k  is  the  local  surface  curvature. 
On  the  gas  side  of  the  interface,  Bernoulli’s  equation  is: 


+  P,  =  0  (10) 

Mathematically,  Eqs.  8-10  provide  a  well-posed  set  of  equations  for  velocity  potentials,  the  gas 
pressure  at  the  interface,  and  the  shape  of  the  interface  (implicit  in  the  surface  curvature,  /c).  This 
set  of  equations  is  solved  using  a  Boundary  Element  Method  (BEM)  which  begins  with  an  integral 
representation  of  Laplace’s  equation: 

j^[4>^-qG]dT  =  0  (11) 


where  0(n)  is  the  value  of  the  potential  at  a  point  F  denotes  the  boundary  of  the  domain,  and  G 
is  the  free-space  Green’s  function  corresponding  to  Laplace’s  equation.  An  analogous  form  of  Eq. 

11  can  also  be  derived  for  the  gas  phase  potential.  For  a  well-posed  problem,  either  or  g  =  d<j>fdn 
must  be  specified  at  each  “node”  on  the  boundary.  Here  n  is  the  outward  normal  to  the  boundary 
so  that  q  represents  the  velocity  normal  to  the  boundary.  The  quantity  a  in  Eq.  11  results  from 
singularities  introduced  as  the  integration  passes  over  the  boundary  point,  ?%•. 

Using  this  methodology,  models  have  been  developed  for  both  two-dimensionaP^  and  axisymmetric^^ 
flowfields.  If  we  let  r  and  z  denote  radial  and  axial  coordinates,  respectively,  and  denote  the  base 
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the  Green  s  function  solution  to  the  axisymmetric  Laplacian  can  be  writ- 


point  with  subscript 
ten: 


where 


4rA;(p) 

\/(r  +  r,)'^  +  (r- 

(r-r,)^  +  (;- 
(r  +  r,)2  +  (--  -)2 


(12) 

(13) 


and  Keip)  is  the  complete  elliptic  integral  of  the  first  kind.  For  computational  efficiency,  this 
quantity  is  calculated  using  a  curve  fit^^  which  has  an  accuracy  to  10“®. 

In  the  case  of  a  2-D  flow  (letting  x  and  y  represent  the  coordinates),  we  have: 


(14) 


In  both  cases,  we  presume  that  both  0  and  q  vary  linearly  along  the  length  of  a  given  element. 
This  assumption  permits  the  construction  of  a  set  of  matrices  involving  the  nodal  values  of  0 
and  q  and  the  integrals  (over  a  given  element)  given  in  Eq.  11.  For  the  2-D  flows,  integration 
across  a  segment  can  be  carried  out  analytically.  Singularities  resulting  from  integration  across 
a  segment  containing  the  base  point  are  also  integrable^®.  In  the  case  of  axisymmetric  flow,  the 
integrations  must  be  carried  out  numerically^^.  In  this  case,  we  choose  a  four-point  Gaussian 
quadrature  for  evaluation  of  integrals.  Logarithmic  singularities  which  arise  in  the  elliptic  integral 
when  the  segment  contains  the  base  point  are  treated  with  a  special  Gaussian  integration  designed 
to  accurately  treat  this  condition. 

Nodes  on  the  interface  are  assumed  to  travel  with  the  local  liquid  surface  velocity,  so  a  trans¬ 
formation  from  the  Eulerian  to  Lagrangian  reference  frame  is  required.  After  this  transformation, 
Eqs.  9  and  10  become: 


D<t> 

Dt 


(15) 


P9  =  -|(V<^,)'  -  +  eV<^  •  V4>g  (16) 

In  these  expressions,  the  notation  D  /Dt  denotes  changes  in  time  for  nodes  moving  with  the  liquid 
interfacial  velocity. 

Surface  slope  and  curvature  are  obtained  from  a  4th-order  treatment  to  insure  accurate  resolu¬ 
tion  of  the  surface.  Cubic  splines  are  used  to  fit  current  locations  of  surface  nodes  for  the  purpose 
of  regridding  as  the  calculation  proceeds,  Regridding  is  necessary  in  many  calculations  because  of 
the  tendency  of  nodes  to  “bunch”  in  regions  of  highest  curvature.  Models  have  been  validated  (for 
nonlinear  calculations)  to  insure  that  solutions  are  insensitive  to  the  grid  which  has  been  selected. 
Droplets  are  assumed  to  be  “pinched”  from  the  main  body  of  fluid  if  a  node  lies  within  5%  of  the 
orifice  radius  of  the  centerline  (or  another  node).  Numerous  validations  have  proven  that  solutions 
are  insensitive  to  this  “pinch  criterion”.  Using  this  methodology,  the  BEM  formulation  permits 
accurate  solution  of  highly-distorted  interfaces  in  an  unsteady  flow. 


Computational  Grids  and  Boundary  Conditions 
Transverse  Mode  Simulations 

Figure  2  highlights  the  computational  grid  and  boundary  conditions  employed  for  coupled,  2~D 
simulations  of  a  liquid  column  subjected  to  an  acoustic  oscillation.  For  this  problem,  the  Weber 
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number  {We),  gas/liquid  density  ratio  (e),  and  gas/liquid  frequency  ratio  {ug/uj)  are  input  param¬ 
eters  for  a  given  simulation.  Assuming  quantities  are  nondimensionalized  using  the  liquid  density, 
column  radius,  and  peak  acoustic  velocity,  the  dimensionless  liquid  column  natural  frequencies  (Eq. 
5)  can  be  written: 


2  _  en{n^  -  1) 
”  {l  +  e)We 


(17) 


We  presume  that  the  acoustic  oscillation  can  be  represented  by  a  simple  sine  wave,  so  that  the 
gas-phase  velocity  potential  far  from  the  column  may  be  written: 


(j>g  =  xsm{ij:gt/2)  (18) 

The  factor  of  1/2  is  included  inside  the  sine  function  to  account  for  the  fact  that  the  inviscid  solution 
is  insensitive  to  the  direction  of  gas  flow.  Numerical  experiments  indicate  that  farfield  conditions 
may  be  accurately  assumed  if  the  outer  gas  boundary  is  placed  15  jet  radii  from  the  origin. 

For  nodes  along  the  symmetry  axis  in  Fig.  2,  the  normal  velocity,  g,  is  assumed  to  be  zero, 
while  the  Bernoulli  conditions  (Eqs.  15  and  16)  provide  the  necessary  boundary  information  for 
liquid  and  gas  nodes  lying  on  the  free  surface  boundary.  A  stable,  accurate,  treatment  of  the 
coupling  of  the  two  flows  through  the  gas  pressure  (which  appears  in  both  Eqs.  15  and  16)  has 
been  developed^^’^®.  The  treatment  begins  by  solving  for  liquid  surface  velocities  from  Eq.  11 
and  noting  that  Qg  =  —q  on  the  surface  since  we  are  tracking  nodes  with  respect  to  the  motion  of 
the  liquid.  Given  the  qg  value  on  the  surface,  solution  of  the  analagous  form  of  Eq.  11  provides 
current  values  for  <f>g  on  the  interface.  The  current  gas  pressure  is  then  computed  from  Eq.  16 
by  approximating  the  D<))g/Dt  term  using  a  first-order  backward  difference.  Finally,  new  values 
of  <t)  on  the  interface  (for  the  next  timestep)  are  obtained  via  integration  of  Eq.  15.  Simulations 
presented  in  the  following  section  employ  17-33  nodes  along  the  interface  and  typically  run  in  a  few 
hours  on  a  Sun  Sparcstation  5  using  a  dimensionless  time  step  of  0.005. 


Longitudinal- Mode  Simulations 

Longitudinal-mode  simulations  are  performed  by  assuming  that  an  oscillation  in  chamber  pressure 
at  the  injector  face  will  lead  to  a  dynamic  massflow  through  the  orifice.  Under  this  assumption, 
we  employ  the  axisymmetric  model  of  Hilbing,  et  al.^^  with  an  unsteady  inflow  cis  indicated  in  Fig. 
3.  For  this  calculation,  the  presence  of  the  gas  is  neglected,  and  nodes  placed  along  the  orifice  exit 
plane  are  held  fixed.  Along  this  inflow  boundary,  the  velocity  is  assumed  to  be: 

9  =  1  +  9'sin(w5<)  (19) 

where  q'  is  the  fractional  change  (in  magnitude)  of  the  velocity  due  to  pressure  oscillations  at  the 
injector  face.  The  oscillation  is  assumed  to  be  at  a  frequency  identical  to  that  of  the  acoustic 
disturbance,  and  velocities  are  nondimensionalized  against  the  mean  velocity  exiting  the  orifice. 

For  nodes  lying  on  the  free  surface,  Eq.  11  (with  Pg  =  0)  is  integrated  in  time  to  give  4)  values 
along  the  interface  which  serve  as  boundary  conditions.  Gravity  is  neglected  in  the  simulations  since 
its  influence  is  very  small  for  high-speed  jets.  Typical  solutions  employ  a  grid  spacing  corresponding 
to  20%  of  the  orifice  radius.  Run  times  can  vary  substantially  (from  a  few  hours,  to  several  days) 
depending  on  the  length  of  the  jet,  i.e.,  long  jets  have  a  large  number  of  nodes  and  require  more 
calculations  per  timestep. 
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Results 

Transverse  Mode  Simulations 

For  all  simulations  presented  here,  a  gas/liquid  density  ratio  of  e  =  0.01  has  been  assumed.  This 
value  is  typical  of  high  pressure  combustion  conditions  as  can  be  noted  in  Table  1.  Parametric 
studies  involving  this  parameter^^  indicate  that  jet  behavior  is  not  terribly  sensitive  to  e  for  low 
density  ratios. 

Under  this  assumption,  the  response  of  a  column  to  an  acoustic  oscillation  at  the  lowest-order 
natural  frequency  of  the  jet  [ug  =  uj\n=2)  is  shown  in  Fig.  4  for  a  W’eber  number  of  0.1.  As 
expected,  we  see  a  strong  response  under  these  conditions,  with  the  distortion  increasing  with  each 
successive  period  of  the  oscillation.  However,  for  long  times  (several  periods),  a  very  interesting 
behavior  results.  This  behavior  is  best  observed  by  plotting  the  time  history  of  the  position  of  a 
node  on  the  top  of  the  column,  as  shown  in  Fig.  5.  Here  we  note  that  the  overall  magnitude  of  the 
distortion  is  bounded  even  though  there  is  no  dissipation  in  this  inviscid  simulation. 

The  explanation  for  this  phenomena  lies  in  the  nonlinear  behavior  of  the  column  itself.  For  finite 
deformations,  the  natural  frequency  of  oscillation  is  reduced  from  that  predicted  by  linear  theory. 
This  effect  is  well  documented  for  oscillations  of  liquid  droplets,  but  has  not  been  investigated 
substantially  for  liquid  columns.  Figure  6  shows  the  nonlinear  frequency  shift  as  a  function  of 
initial  surface  deflection  for  both  droplets  and  columns.  As  the  jet  distorts  due  to  excitation  at  its 
linear  natural  frequency,  further  excitation  at  this  frequency  will  actually  become  destructive  for 
some  level  of  column  deformation.  Beyond  this  point,  the  oscillation  will  decay  (as  shown  in  Fig. 
5)  until  the  cycle  repeats. 

The  overall  frequency  response  of  the  column  is  shown  in  Fig.  7.  As  expected,  we  see  a  strong 
peak  at  excitation  frequencies  near  the  natural  harmonic  (n  =  2)  frequency  of  the  column.  However, 
it  is  quite  interesting  to  note  that  very  little  response  is  obtained  for  all  other  frequencies  in  the 
domain.  There  is  a  small  peak  at  the  subharmonic  Ug  =  but  there  is  very  little  activity  at 

the  higher  order  (4th,  6th,  etc.).  We  expect  very  little  response  ifcj^  >>  u  since  the  jet  inertia  does 
not  permit  reactions  on  this  short  time  scale.  In  the  other  limit  Ug  «  uj,  the  column  response  is 
nearly  quasi-steady  and  peak  deflections  are  simply  a  function  of  We.  Note  that  at  this  condition 
(We  =  0.1)  very  little  jet  broadening  would  be  predicted  for  a  steady  crossflow  cls  indicated  in  the 
far  left  half  of  the  curves  in  Fig.  7. 

These  results  have  obvious  implications  to  LRE  combustion  instability.  A  design  which  has 
transverse  acoustic  modes  near  the  jet  natural  frequency  (n=2  in  Eq.  6)  will  be  subject  to  coupling 
from  wave  structures  present  in  the  gets  domain.  If  the  jet  broadens  substantially,  its  drag  will 
increase  and  it  will  be  more  subject  to  deflection  by  the  transverse  flow.  For  impinging  element  de¬ 
signs,  the  broadening  will  have  direct  consequences  on  the  atomization/mixing  in  the  impingement 
region.  As  indicated  in  Table  1,  the  F-1  injector  natural  frequencies  do  lie  near  the  range  of  tan¬ 
gential  acoustic  frequencies  experienced  during  testing  of  this  engine.  Therefore,  this  phenomenom 
could  explain  (at  least  in  part)  the  tangential  stability  problems  encountered  by  this  engine.  The 
“good  news”  is  that  this  phenomenom  appears  to  be  limited  to  the  lowest-order  natural  frequency 
of  the  column;  subharmonics  and  higher  harmonics  do  not  play  a  substantial  role  as  indicated  in 
Fig.  7. 

Longitudinal  Mode  Simulations 

Currently,  simulation  of  a  three-dimensional,  viscous,  turbulent,  two-phase  flow  which  corresponds 
to  “real  world”  conditions  is  well  beyond  current  modeling  capabilities.  For  this  reason,  the  longi¬ 
tudinal  simulations  have  been  limited  to  a  relatively  low  velocity  regime  in  which  the  assumptions 
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of  axisymmetry  and  inviscid  flow  are  most  appropriate.  The  main  focus  of  these  simulations  is  to 
address  the  effects  of  unsteady  perturbations  on  the  character  of  breakup  under  these  assumptions. 
Since  the  simulations  in  this  section  do  not  include  the  influence  of  gas-phase  pressure  variations 
about  the  jet,  the  Weber  number  is  defined  in  terms  of  the  liquid  density.  A  Weber  number  of  100 
is  assumed  for  results  presented  herein. 

Figures  8  and  9  address  the  effect  of  the  effect  of  disturbance  magnitude  on  the  atomization 
process  for  a  disturbance  at  the  jet  natural  frequency.  In  this  case,  analysis  of  Eq.  6  indicates  that 
a  dimensionless  frequency  of  u;  0.7  generates  waves  of  the  most  unstable  length.  In  Fig.  8,  one 
can  note  that  nonlinear  effects  lead  to  the  formation  of  a  “main”  and  “satellite”  droplet  from  a 
single  wavelength  of  the  instability.  This  behavior  is  well  known  for  low-speed  jets.  Periodic  bulges 
appear  in  the  jet  due  to  the  unsteady  inflow,  as  noted  experimentally  by  Reba  and  Brosilow^. 

Increctsing  the  size  of  the  perturbation  not  only  decreases  jet  breakup  length,  but  also  effects 
the  shape  and  size  of  both  main  and  satellite  drops.  As  g'  is  increased,  main  droplets  take  on  a 
“squashed”  shape  as  a  result  of  high  velocity  fluid  (when  sin(a;^)  >  0)  encountering  lower  velocity 
fluid  which  has  already  exited  the  nozzle.  This  phenomena,  known  as  the  Klystron  effect,  has  been 
documented  qualitatively  by  numerous  researchers^ Quantitative  comparisons  of  droplet  sizes 
are  compared  with  an  inflnitessimal  disturbance^^  in  Fig.  9,  Note  that  the  inviscid  predictions 
match  experiments  quite  well  in  this  regime.  Viscosity  tends  to  slow  down  the  breakup  process, 
but  does  not  fundamentally  change  droplet  sizes  in  this  flow  regime.  The  satellite  size  is  shown 
to  increase  with  perturbation  amplitude  at  the  frequency  ratio  selected  for  these  simulations.  Ex¬ 
trapolation  to  higher  q'  values  would  presumably  lead  to  a  monodisperse  case  in  which  both  drops 
are  the  same  size.  Through  the  use  of  piezoelectric  drivers,  droplet  trains  of  this  type  have  been 
created  experimentally^^. 

Figures  10  and  11  address  the  effect  of  disturbance  frequency  on  the  character  of  the  jet.  As  in 
the  case  of  amplitude  dependence  (Figs.  8,9),  the  jet  behavior  is  also  frequency  dependent.  In  this 
case,  increasing  disturbance  frequency  tends  to  decrease  breakup  length  even  beyond  >  1; 

a  trend  not  predicted  by  linear  theory.  In  addition,  the  size  of  satellite  drops  tends  to  decrease 
with  increasing  frequency  as  shown  in  Fig.  10.  At  a  frequency  ratio  of  1.57,  satellite  drops  have 
nearly  vanished  indicating  that  a  frequency  near  this  value  can  produce  monodisperse  atomization. 
In  fact,  viscous  effects  would  probably  preclude  the  formation  of  this  very  small  structure.  Quan¬ 
titative  predictions  of  droplet  sizes  are  shown  in  Fig.  11.  Here,  the  solid  lines  reflect  results  for 
an  infinitessimal  perturbation,  as  in  Fig.  9  and  the  datapoints  are  measurements  of  Rutland  and 
Jameson^^  and  Lafrance^^. 

A  final  simulation  was  conducted  for  the  case  of  a  high  Weber  number  {We  =  1000),  large 
amplitude  (g'  =  10%),  perturbation.  A  time  sequence  for  this  case  is  shown  in  Fig.  12.  In  this 
case,  the  Klystron  effect  leads  to  a  large  radial  broadening  of  the  jet,  cts  observed  in  numerous 
experiments^ Large  perturbations  of  this  nature  have  obvious  repercussions  to  jet  behavior, 
causing  clustered  “clumps”  of  fluid  at  locations  corresponding  to  peak  flowrates  during  the  unsteady 
process.  In  addition,  the  droplets  formed  from  the  annular  rings  of  fluid  shown  in  Fig.  12  will  be 
much  smaller  than  drops  formed  from  a  “steady”  atomization  process;  an  effect  observed  by  Ingebo^. 
Additional  modeling  efforts  will  be  required  to  obtain  a  more  quantitative  evaluation  of  this  effect. 


Conclusions 

Nonlinear  two-dimensional  and  axisymmetric  numerical  simulations  have  been  applied  to  simulate 
the  effect  of  acoustic  perturbations  on  liquid  jet  atomization  processes.  A  liquid  column  is  shown  to 
react  very  strongly  to  perturbations  at  its  lowest-order  natural  frequency  uj  =  \/(6cr)/((/?+  Pg)(i^), 


41 


but  has  minimal  unsteady  response  for  frequencies  not  near  this  value.  Early  (unstable)  fuel 
injector  designs  used  in  the  F-1  engine  program  were  shown  to  be  near  this  natural  frequency.  In 
addition,  nonlinear  natural  frequencies  for  columns  under  finite  deformation  have  been  quantified. 
The  frequency  shift  associated  with  a  finite  amplitude  deformation  is  shown  to  cause  a  limiting 
amplitude  of  column  distortion. 

The  response  of  a  liquid  jet  to  longitudinal  excitation  (as  modeled  through  the  use  of  a  dynamic 
orifice  massflow)  has  also  been  quantified.  Periodic  bulges  in  the  jet’s  surface,  typically  referred  to  as 
the  ''Klystron  effect”,  are  described  quantitatively  using  this  model.  For  low  speed  jets,  increasing 
both  amplitude  and  frequency  of  the  disturbance  is  shown  to  increase  the  size  of  "satellite”  drops 
formed  by  nonlinear  deformation  of  the  column.  At  large  excitation  amplitude,  distinct  "mushroom¬ 
shaped”  structures  appear,  as  described  by  several  researchers.  It  is  obvious  that  the  formation  of 
these  structures  will  decrease  mean  drop  size,  but  more  effort  is  required  to  quantify  the  extent  of 
the  size  change. 
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Figure  Captions 


1.  Liquid  Jets  Subjected  to  Transverse  and  Longitudinal  Acoustic  Oscillations 

2.  Computational  Domain  and  Boundary  Conditions  for  Transverse  Mode  Simulations 

3.  Computational  Domain  and  Boundary  Conditions  for  Longitudinal  Mode  Simulations 

4.  Column  Shapes  at  Various  Times  of  an  Acoustic  Perturbation,  Ug  =  lj.  We  =  0.1,  e  = 

0.01 

5.  Motion  of  a  Node  on  the  Top  of  the  Column  During  Acoustic  Perturbation,  Wj  =w.  We  = 

0.1,  e  =  0.01 

6.  Nonlinear  Frequency  Shift  for  Liquid  Drops  and  Columns  under  Vaccuum  (or  Low  Gas 
Density)  Conditions 

7.  Nonlinear  Frequency  Response  of  Liquid  Column;  Maximum  (Prolate)  and  Minimum 
(Oblate)  Aspect  Ratios  Obtained  During  Oscillation 

8.  Effect  of  Longitudinal  Disturbance  Amplitude  on  Behavior  of  Liquid  Jet  at  IVe  =  100,  = 

w  =  0.7.  (a)  q'  =  2%,  (b)  q'  =  4%,  (c)  q'  =  6% 

9.  Effect  of  Longitudinal  Disturbance  Amplitude  on  Drop  Size  for  We  —  100,  Ug  =  u)  =  0.7. 
Solid  Lines  are  for  Infinitessimal  {q'  «  0)  Disturbance. 

10.  Effect  of  Longitudinal  Disturbance  Frequency  on  Behavior  of  Liquid  Jet  at  We  =  100,  q'  = 
2%.  (a)  tj/ujg  =  0.71,  (b)  uj/u)g  =  1.0,  (c)  uj/u)g  =  1.29,  (d)  u)/ug  —  1.57 

11.  Effect  of  Longitudinal  Disturbance  Frequency  on  Drop  Size  for  Liquid  Jet  at  We  = 
100,  q'  =  2%.  Solid  Lines  are  for  q'  v  0. 

12.  Jet  Behavior  Under  Violent  Oscillation,  We  =  1000,  ^  =  1,  q'  =10%. 
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Figure  1:  Liquid  Jets  Subjected  to  Transverse  and  Longitudinal  Acoustic  Oscillations 
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Figure  2;  Computational  Domain  and  Boundary  Conditions  for  Transverse  Mode  Simulations 
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Figure  4:  Column  Shapes  at  Various  Times  of  an  Acoustic  Perturbation,  Ug  =  ui,  lVe  =  0.1,  e  = 
0.01 
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Figure  7:  Nonlinear  Frequency  Response  of  Liquid  Column;  Maximum  (Prolate)  and  Minimum 
(Oblate)  Aspect  Ratios  Obtained  During  Oscillation 
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Figure  8:  Effect  of  Longitudinal  Disturbance  Amplitude  on  Behavior  of  Liquid  Jet  at  We 
100,  Wg  =  w  =  0.7.  (a)  q'  =  2%,  (b)  g'  =  4%,  (c)  q'  =  6% 
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Figure  9:  Effect  of  Longitudinal  Disturbance  Amplitude  on  Drop  Size  for  We 
Solid  Lines  are  for  Infinitessimal  (q'  ss  0)  Disturbance. 


Figure  10:  Effect  of  Longitudinal  Disturbance  Frequency  on  Behavior  of  Liquid  Jet  at  We 
100,  q'  =  2%.  (a)  uj/ujg  —  0.71,  (b)  uj/uig  =  1.0,  (c)  =  1.29,  (d)  uifwg  =  1.57 


9  Appendix  D  -  Acoustic  Interaction  with  a  Droplet 
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Abstract-The  unsteady,  nonlinear  response  of  a  liquid  droplet  to  an  imposed  acoustic  perturbation  has 
been  simulated  using  a  model  based  on  Boundary  Element  Methods  (BEMs).  The  model  was  used  to  study 
the  influence  of  the  acoustic  frequency,  intensity,  and  gas/liquid  density  ratio  upon  the  droplet  behavior. 
Heightened  droplet  responses  were  observed  for  frequencies  near  the  harmonics  of  the  second  and  fourth  mode 
frequencies.  Several  types  of  droplet  atomization  have  been  observed  as  the  acoustic  intensity  is  increased. 
Increasing  gas  density  (at  fixed  excitation  conditions)  also  heightens  droplet  response. 

Key  Words:  drop  break  up,  atomization,  drop  dynamics,  acoustics,  drop  oscillations 


1  Introduction 

The  behavior  of  liquid  droplets  in  the  presence  of  an  acoustic  field  is  a  phenomena  of  fundamental  importance. 
There  are  a  variety  of  applications  in  which  droplets  are  excited  by  acoustic  energy,  such  as  meteorological 
physics,  containerless  processing,  and  atomization.  Droplet  oscillations  in  the  absence  of  external  excitation 
will  ultimately  dissipate,  and  the  drop  will  return  to  its  equilibrium  state.  However,  a  strong  acoustic  field 
may  be  present,  such  is  the  case  in  airborne  combustors  and  acoustic  levitators.  In  such  a  situation,  the 
droplet  will  experience  forced  oscillations.  This  oscillatory  behavior  of  drops  may  greatly  affect  the  process 
of  secondary  atomization. 

A  large  number  of  researchers  have  studied  “free  oscillations”  of  droplets  under  conditions  where  forces 
from  the  gaseous  phase  are  neglected.  The  initial  linear  analyses  of  this  problem  are  due  to  Rayleigh  (1879) 
and  Lamb  (1932)  for  the  inviscid  and  weak- viscous  cases.  These  results  provide  the  droplet’s  frequency  of 
oscillation  under  various  modes  under  the  presumption  that  the  shape  perturbation  is  infinitessimal.  Lamb’s 
result  indicates  that  weak  viscous  effects  have  a  negligible  effect  on  the  frequency  predicted  by  Rayleigh. 
Prosperetti  (1980)  furthered  the  linear  analysis  for  arbitrary  viscous  effects.  Despite  his  improvements  to 
the  linear  analysis,  there  was  considerable  disparity  between  his  predicted  decay  rate  of  oscillations  and  the 
experimentally  observed  rate. 

More  recently,  theoretical/numerical  efforts  have  focused  on  nonlinear  effects.  Tsamopoulos  and  Brown 
(1983)  developed  a  theoretical  series  solution  for  moderate  amplitude,  inviscid  droplet  oscillations.  They 
determined  that  the  droplet  resonant  frequency  decreases  with  the  square  of  the  oscillation  amplitude. 
Lundgren  and  Mansour  (1988)  used  a  boundary  integral  method  to  develop  a  model  for  large  amplitude 
oscillations  of  droplets  in  zero-gravity  including  weak  viscous  effects.  They  discovered  that  relatively  small 
viscosities  can  significantly  affect  the  coupling  of  oscillatory  modes.  However,  their  method  was  limited  to 
only  weak  viscous  problems.  More  recently,  Basaran  (1991)  used  the  Galerkin/finite  element  technique  to 
model  nonlinear  oscillations  of  viscous  droplets.  From  these  simulations  he  observed  that  a  finite  viscosity 
has  a  much  larger  effect  on  mode  coupling  than  what  is  predicted  by  the  calculations  including  weak  viscous 
effects. 

Some  notable  experimental  work  in  droplet  oscillations  has  been  conducted  by  Trinh  and  Wang  (1981), 
in  which  large  amplitude,  nonlinear  oscillations  of  drops  were  studied.  In  their  work,  a  neutrally-buoyant 
droplet  was  suspended  in  an  immiscible  liquid  and  excited  by  acoustic-radiation  pressure  forces  generated 
by  an  acoustic  levitator.  More  recently,  Wang,  Anilkumar,  and  Lee  (1996)  studied  the  oscillations  of  low- 
viscosity  drops  in  a  microgravity  environment  on  board  the  space  shuttle  using  a  similar  acoustic  chamber 
to  induce  droplet  deformation.  In  these  experiments,  droplets  were  deformed  and  then  allowed  to  oscillate 
freely.  From  a  controlled  break  up  of  a  liquid  jet,  Becker,  Hiller,  and  Kowalewski  (1990)  generated  virtually 
monodispersed  droplets  that  oscillated  in  a  damped,  axisymmetric  fashion.  These  experiments  verified  the 
reduction  in  natural  frequency  at  finite  amplitudes  as  predicted  by  Tsamopolous  and  Brown. 

Despite  these  advancements  in  the  area  of  free-oscillations  of  droplets,  there  have  been  relatively  few  works 
aimed  at  increasing  the  understanding  of  droplet  response  to  forced  oscillations.  Here,  the  experimental  work 
of  Daidzic  (1995)  is  a  notable  exception.  Using  an  acoustic  levitator,  Daidzic  examined  nonlinear  forced 
oscillations  of  droplets.  In  his  experiments  the  droplets  exhibited  “chaotic  behavior” ,  presumably  because 
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the  forcing  function  was  three-dimensional  and  time-dependent.  He  concluded  that  prediction  of  droplet 
behavior  over  a  long  period  of  time  is  difficult  at  best,  and  warrants  further  investigation. 

The  works  of  the  aforementioned  researchers  have  considerably  improved  the  understanding  of  droplet 
behavior.  However,  the  forced  oscillation  problem  still  warrants  additional  consideration.  The  focus  of 
this  work  is  to  examine  forced,  nonlinear,  inviscid  droplet  oscillations  by  computational  analysis.  The 
Boundary  Element  Method  (BEM)  will  be  used  to  conduct  these  studies.  This  technique  offers  a  savings  in 
computational  power  as  compared  to  other  methods,  such  as  computational  fluid  dynamics,  while  maintaining 
a  high  degree  of  accuracy.  As  the  name  implies,  BEMs  require  the  discretization  of  only  the  boundary  of  the 
domain,  providing  a  drastic  reduction  in  the  total  number  of  nodes  (as  compared  to  a  mesh  based  schemes) 
needed  for  an  accurate  solution.  Hilbing,  Heister,  and  Spangler  (1995)  and  Mansour  and  Lundgren  (1990) 
have  also  demonstrated  a  capability  of  running  BEM  simulations  beyond  atomization  events.  The  following 
sections  of  this  paper  describe  the  development,  validation,  and  results  generated  using  this  computational 
tool. 


2  Model  Development 


Under  many  situations,  the  wavelength  of  the  acoustic  perturbation  is  much  greater  than  the  droplet  radius, 
which  implies  that  spatial  variations  within  the  acoustic  wave  are  negligible  and  that  the  disturbance  can 
be  modeled  as  an  unsteady,  incompressible  flow.  Further,  we  assume  an  axisymmetric  domain  and  neglect 
viscosity  in  both  gas  and  liquid  phases.  Under  these  conditions,  the  dynamics  of  both  liquid  and  gas  phases 
are  described  by  Laplace’s  equation: 


V“(/>  =  =  0 


(1) 


where  (j>  and  (i)g  are  velocity  potentials  for  liquid  and  gaseous  phases,  respectively. 

If  we  choose  the  droplet  radius  (a),  peak  speed  in  the  acoustic  disturbance  [U],  and  liquid  density  [p] 
as  dimensions,  the  interaction  between  the  droplet  and  the  acoustic  disturbance  is  characterized  by  the 
gas/liquid  density  ratio, 


(2) 


the  Weber  number  batsed  on  gas  density, 


(3) 


and  the  frequency  ratio, 


w 


(4) 


Changes  in  the  magnitude  of  the  acoustic  disturbance  are  introduced  through  the  Weber  number,  which  is 
the  ratio  of  the  aerodynamic  forces  to  surface  tension  (cr)  forces.  The  frequency  ratio  is  the  ratio  of  the 
acoustic  frequency  of  the  gas,  a;,  to  the  linear  natural  frequency  of  the  second  mode  for  a  liquid  drop,  In 
the  following  development,  we  presume  that  the  nondimensionalization  described  above  has  been  applied. 

The  droplet’s  fundamental  frequencies  are  obtained  from  the  classic  analysis  by  Lamb.  The  nondimen- 
sional  form  for  the  frequency  of  mode  “m”  is: 


m(m  -I-  l)(m  —  l)(m  -h  2) 
(m  -h  1  +  me)  We 


(5) 


For  this  case,  the  lowest-order  (m  —  2),  or  natural  frequency  of  a  droplet  reduces  to: 

2  _ 

"2  =  (3  +  2£)We 


(6) 


We  expect  strong  droplet  response  when  the  acoustic  excitation  frequency  uj  lies  near  harmonics  of  u/n- 
At  any  instant  in  time,  Eqn.  1  provides  a  connection  between  values  of  the  velocity  potentials  and 
velocities  measured  normal  to  the  local  surface  [q  and  qg)  for  liquid  and  gaseous  phases,  respectively.  This 
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equation  is  solved  using  the  BEM  by  beginning  with  its  integral  representation.  For  the  liquid  phase,  the 
resulting  integral  representation  of  Eqn.  I  becomes: 


(/r  =  o 


(7) 


where  0[f])  is  the  velocity  potential  at  a  point  n,  F  is  the  boundary  of  the  domain,  and  G  denotes  the 
free  space  Green's  function  corresponding  to  Laplace's  equation.  In  addition,  h  is  taken  to  be  the  outward 
normal  to  the  domain  boundary,  and  a  results  from  singularities  introduced  as  the  integration  passes  over 
the  boundary  point,  r,.  By  using  Eqn.  7,  it  is  possible  to  solve  for  (p  ov  q  provided  that  one  of  them  is 
known  at  each  point  on  the  boundary.  If  we  let  r  and  c  represent  radial  and  axial  coordinates  and  use  the 
subscript  i  to  denote  the  “base  point”  where  the  integration  takes  place,  the  free  space  Green’s  function  for 
the  axisymmetric  Laplacian  can  be  expressed  (Liggett  and  Liu  1983): 


G  = 


4rK{p) 

Vir+ nr- +  (=-=<)- 


(8) 


where 


P  “ 


(r-r,)^  +  (z-r,)^ 

(r  +  r,)2  +  (r-;,)2 


(9) 


and,  K(p)  is  the  complete  elliptic  integral  of  the  first  kind.  For  computational  efficiency,  this  parameter  is 
curvefit  (to  an  accuracy  of  10“®)  using  results  from  Abramowitz  and  Stegun  (1970). 

The  integration  in  Eqn.  7  is  performed  by  discretizing  the  boundary  into  a  finite  number  of  segments. 
Along  each  segment,  both  (f)  and  q  are  assumed  to  vary  linearly.  Integrations  are  carried  out  by  letting  each 
node  on  the  boundary  represent  a  base  point,  yielding  a  set  of  linear  equations  relating  4>  and  q  involving 
all  boundary  nodes.  The  fully**populated  matrices  used  to  store  coefficients  of  these  equations  are  inverted 
using  the  Grout  Method  LU  Decomposition  from  Numerical  Recipes  in  Fortran  (Beyer  1991).  Additional 
details  regarding  the  BEM  solution  procedure  are  provided  in  Hilbing,  Heister,  and  Spangler  (1995). 


2.1  Free  Surface  Treatment 


A  procedure  similar  to  that  of  Longuet-Higgins  and  Cokelet  (1976)  is  used  to  update  the  position  of  nodes 
on  the  interface.  Free  surface  nodes  are  “tracked”  along  lines  parallel  with  the  local  velocity  vector  in  the 
liquid.  Under  this  assumption,  flow  kinematics  require: 


Dr 

'm 


dr 


Dz 

Dt 


dz 


(10) 


The  velocities  calculated  by  solving  Laplace’s  equation  are  normal  and  tangential  to  the  surface.  While 
the  normal  velocity  is  generated  via  the  solution  of  Laplace’s  equation,  the  tangential  velocity  {d<j)/ds)  is 
calculated  using  a  4th-order  (five  point)  centered  difference  method.  Radial  and  axial  velocities  required  in 
Eqn.  10  are  determined  via  the  coordinate  transformation: 


and 


d<i>  d<f>  .  .  . 

—  =  —  sin/? +  ^  cos/? 


dcf>  d<t>  .  .  . 

_  =  _cos/?-gsm/? 


(11) 


(12) 


where  /?  is  the  local  wave  slope. 

The  unsteady  Bernoulli  equation  provides  the  dynamic  boundary  condition  for  nodes  on  the  interface.  In 
an  Eulerian  system  where  time  derivatives  are  assumed  to  occur  at  a  fixed  spatial  location,  the  dimensionless 
form  of  this  relation  valid  in  the  liquid  domain  is: 


dt 


+  5(V«'  +  P,  +  ^=0 


and  the  gas  domain  analog  is: 


(13) 

(14) 
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Since  the  nodes  on  the  interface  are  assumed  to  travel  with  the  local  liquid  surface  velocity,  a  transfor¬ 
mation  from  the  Eulerian  to  Lagrangian  reference  frame  is  required.  Under  this  cissumption,  the  following 
equation  can  be  used  to  perform  the  Eulerian  -  Lagrangian  transformation  [23]: 


D(  ) 

Dt 


dt 


+  V0-V{) 


By  applying  this  transformation,  the  Bernoulli’s  equations  (13,14)  become; 


(15) 


D6 

Dt 


2(Vd)- 


(16) 


Dt 


+  6  V0  •  V  <j)g 


(17) 


where  D{  )/Dt  denotes  changes  in  time  for  nodes  moving  with  the  liquid  interface  velocity. 

In  the  case  of  the  liquid  domain,  the  substantial  derivative,  D(t>/ Dt,  is  calculated  using  Eqn.  16.  A  first 
order  backward  difference  method  was  used  to  calculate  the  gas  based  substantial  derivative  appearing  in 
Eqn.  17: 


Dt  At 


(18) 


where  n  indicates  the  time  level. 

A  stable,  time-accurate  procedure  has  been  developed  to  advance  the  solutions  of  (16,  17)  for  all  free- 
surface  nodes.  The  following  steps  are  taken: 

1.  An  initial  value  of  0  and  (pg  are  given  at  each  node  on  the  interface. 

2.  Solution  of  Laplace’s  equation  (described  in  previous  section)  provides  the  liquid  domain  velocities,  q, 
at  each  node  on  the  interface. 


3.  Set  qg  =  —q  since  the  gas  nodes  are  fixed  to  move  with  the  liquid  nodes  and  the  outward  normals  are 
in  opposite  directions. 


4.  Solve  Laplace’s  equations  for  the  gas  domain  to  calculate  0^  at  each  node  in  the  gets. 

5.  Calculate  D(j)g/Dt  using  Eqn.  18. 

6.  Calculate  the  gas  pressure  distribution  along  the  interface  using  Eqn.  17. 

7.  This  value  of  gas  pressure  is  used  to  solve  for  0  at  the  new  time  step  via  integration  of  Eqn  16. 

8.  The  interface  is  “regridded”  using  cubic  splines  of  r,  2r,  Pg,  and  0  to  preserve  the  even  spacing  between 
nodes  (Hilbing  et  al.  1995). 

9.  Steps  2-4  are  repeated  with  the  “regridded”  properties,  in  order  to  determine  the  boundary  conditions 
at  the  next  time  step  for  the  new  grid. 


A  4th-order  Runge-Kutta  scheme  is  employed  in  the  time  integrations.  Using  this  scheme,  it  is  necessary 
to  solve  the  Laplace  equation  eight  times,  four  for  the  liquid  and  four  for  the  gas.  However,  at  the  end  of 
each  time  step  it  is  necessary  to  solve  the  Laplace  equation  for  both  domains  again  (due  to  the  regridding), 
therefore  Laplace’s  equation  is  solved  a  total  of  ten  times  per  time  step  using  this  procedure. 


2.2  Domain  Discretization  and  Boundary  Conditions 

The  computational  domain  and  the  boundary  conditions  used  in  this  analysis  are  displayed  in  Figure  1. 
In  this  figure,  0  and  <j)g  represent  the  velocity  potentials  in  the  liquid  and  gas  domains  respectively,  and 
q  —  d(j>fdh  stands  for  the  velocity  in  the  normal  direction.  The  radial  and  axial  directions  are  denoted  by 
r  and  z.  The  flow  is  considered  axisymmetric  with  respect  to  the  z-direction,  so  that  only  a  half  section  is 
needed  for  the  model.  However,  another  axis  of  symmetry  exists  along  the  radial  axis  due  to  the  assumption 
of  inviscid,  incompressible  flow.  As  a  result,  only  a  quarter  section  is  required  for  this  analysis.  Nodes  have 
been  placed  on  the  interface  between  the  domains,  along  the  radial  axis  in  both  domains,  and  the  outer 
boundary  of  the  gas.  The  gas  nodes  are  denoted  with  a  “x”,  while  the  liquid  nodes  are  labeled  with  an  “o”. 
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Along  the  radial  line  of  symmetry  for  both  domains,  7  —  0  in  the  liquid  and  d^gjdv  =  0  in  the  gas. 
Equations  16  and  17  were  used  to  set  the  boundary  condition  along  the  interface.  Along  the  outer  boundary 
of  the  gas  domain,  a  sinusoidal  velocity  history  is  used  to  model  the  acoustic  disturbance: 

0g-  ZgZO^\yji)  (19) 

Boundary  conditions  are  not  specified  along  the  axis  of  axisymmetry,  since  the  Green  s  functions  vanish  as 
r  0.  For  this  reason,  nodes  are  not  required  along  this  boundary. 


3  Model  Implementation  and  Validation 

Once  the  model  had  been  fully  developed,  its  accuracy  was  checked  by  comparing  its  results  to  results  from 
analytic  and  other  numerical  methods.  Three  separate  validation  cases  were  considered:  steady  flow  over  a 
sphere,  non-linear  oscillations  of  a  liquid  droplet,  and  droplet  profiles  for  a  steady  crossflow. 

The  accuracy  of  the  gas  domain  solution  was  tested  by  comparing  the  distribution  of  the  velocity  potential 
over  a  sphere  with  the  analytic  solution  of  the  same  problem.  The  analytic  solution  was  compared  to  the 
numerical  solutions  for  grids  using,  45,  35,  25,  and  15  nodes  along  the  interface.  The  numerical  results 
agree  well  with  the  analytic  solution,  and  the  magnitude  of  the  error  is  well  below  0.5%  for  all  cases.  For 
computational  efficiency,  it  was  decided  to  employ  only  15  nodes  to  model  the  surface  of  the  droplet.  For 
cases  when  the  drop  becomes  highly  deformed,  as  many  as  45  nodes  were  employed  in  order  adequately 
resolve  the  surface. 

As  explained  previously,  an  expression  for  the  fundamental  oscillatory  mode  was  developed  by  Lamb 
using  linear  analysis.  This  analysis  assumed  a  linearized  surface  shape  of  the  form: 

r  =  1  +  7;cos(n^)  sin(u;nO  (20) 

where  7?  is  the  measure  of  the  initial  disturbance  amplitude.  Tsamopoulos  and  Brown  (1983)  later  expanded 
upon  this  analysis  by  including  second  order  effects.  They  observed  that  the  frequency  of  oscillation  decayed 
with  the  square  of  the  initial  deformation.  The  accuracy  of  the  liquid  domain  solution  was  determined 
by  making  comparisons  between  the  oscillatory  frequencies  from  the  numerical  model  and  the  results  of 
Tsamopoulos  and  Brown.  The  numerical  model  was  run  for  an  array  of  initial  surface  distortions  with  23 
nodes  along  the  boundary.  Figure  2  shows  the  comparison  of  the  results  from  the  two  models.  Excellent 
agreement  exists  between  the  two  results.  It  was  concluded  that  the  liquid  domain  solution  is  accurate. 

The  last  section  of  the  validation  process  was  to  check  the  accuracy  of  the  model  when  there  is  a  strong 
coupling  effect  between  the  gas  and  liquid  domains.  This  was  accomplished  by  comparing  droplet  profiles 
calculated  by  the  numerical  model  and  the  profiles  predicted  by  Miksis,  Vanden-Broeck,  and  Keller  (1981). 
They  developed  a  model  to  predict  the  equilibrium  shape  of  a  droplet  subjected  to  a  uniform  flow.  The 
equilibrium  drop  profile  was  determined  by  balancing  the  pressure  of  the  gas  on  the  surface  with  the  surface 
tension  and  internal  pressure  at  the  same  point. 

In  order  to  achieve  a  steady  state  solution  with  the  model,  it  weis  necessary  to  introduce  a  dissipative 
mechanism,  which  ultimately  yielded  a  steady  state  solution.  Numerical  smoothing  was  applied  to  0  until 
surface  velocities  vanished.  Figure  3  shows  a  comparison  between  the  predicted  profiles  of  the  numerical 
model  and  the  method  used  by  Miksis  et  al.  These  comparisons  clearly  show  excellent  agreement  between 
the  two  models.  Hence,  it  was  concluded  that  the  coupled  solution  is  accurate. 

4  Results 


The  model  was  utilized  to  assess  the  influence  of  perturbation  frequency  (w/cj^),  Weber  number  (We), 
and  density  ratio  (e)  on  the  nonlinear  response  of  a  droplet  to  an  imposed  acoustic  oscillation.  The  results 
presented  here  represent  over  200  runs  of  the  model.  Calculations  typically  used  15  nodes  along  the  interface, 
15  along  the  radial  axis  in  the  gas,  and  10  nodes  along  the  outer  boundary  of  the  gas  domain.  However, 
larger  grids  (as  many  as  45  nodes  along  the  interface)  were  utilized  to  resolve  more  complex  droplet  shapes, 
such  as  those  in  Figure  14  described  below. 

During  these  simulations,  the  level  of  the  droplet  response  is  characterized  by  two  parameters,  the  aspect 
ratio  and  the  oscillatory  mode  coefficients.  The  aspect  ratio  provides  a  gross  meatsure  of  overall  droplet 
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deformation  and  is  defined  as  being  the  ratio  of  the  major  axis  over  the  minor  axis  of  the  droplet  at  peak 
deformation.  The  mode  coefficients,  are  determined  by  expressing  the  current  droplet  shape  as  a  sum  of 
Legendre  Polynomials: 

CO 

r  —  I  +  ^  Cm(t)  cos{7n0)  (21) 

m=l 

where  0  is  the  angle  measured  from  the  positive  r  axis  and  Cn  are  the  time  varying  coefficients.  Because 
Cosine  functions  are  orthogonal,  it  is  possible  to  calculate  the  coefficients  using  (Lundgren  and  Mansour 
1988); 

Cm  =  ^  f  ^  I)cos(m0)d9  (22) 

^  J-TT 

Due  to  the  fact  that  both  incompressible  and  inviscid  assumptions  are  employed  in  the  gcis  phase,  there 
is  no  mechanism  to  generate  a  surface  pressure  distribution  which  is  asymmetric  about  the  r  axis  in  Fig.  1. 
For  this  reason,  the  model  cannot  excite  odd  modes  of  oscillation.  This  fact  was  verified  by  demonstrating 
that  mode  coefficients  for  the  odd  modes  were  all  zero  to  within  machine  accuracy  (10“^"^). 

4.1  Frequency  Response  Spectrum 

The  effect  of  the  acoustic  perturbation  frequency  was  investigated  by  performing  approximately  100  simu¬ 
lations.  While  the  frequency  was  varied,  the  Weber  number  was  held  at  0.5779  and  the  density  ratio  was 
0.00123.  These  conditions  correspond  to  acoustic  excitation  of  a  100  micron  water  droplet  with  a  160  db 
disturbance  in  ambient  air.  A  time  step  of  0.05  seconds  was  employed  for  these  computations.  Results  are 
displayed  on  Figure  4,  which  charts  the  aspect  ratio  for  a  range  of  frequency  ratios.  Readers  are  cautioned 
that  the  purpose  of  this  investigation  is  to  determine  regions  of  high  response;  actual  aspect  ratios  do  depend 
on  the  assumed  initial  conditions  (spherical  vs.  deformed  drop,  sine  wave  vs.  cosine  wave  disturbance). 

Figure  4  displays  a  series  of  peaks  that  occur  near  the  harmonics  of  the  natural  frequency  of  the  droplet. 
A  noteworthy  area  is  the  break  up  region  that  exists  for  frequency  ratios  between  0.80  and  0.90.  In  this  band, 
the  disturbance  is  tuned  to  the  natural  frequency  of  the  droplet  to  an  extent  such  that  atomization  occurs. 
All  peaks  are  shifted  to  frequencies  slightly  less  than  the  linear  result  due  to  the  nonlinear  frequency  shift 
as  shown  in  Figure  2.  The  continual  shifting  of  droplet  natural  frequency  with  deformation  level  leads  to  a 
bounded  response  over  much  of  the  frequency  range  (for  this  Weber  number).  In  these  regions,  the  acoustic 
perturbation  continues  to  excite  the  droplet  until  the  oscillations  of  the  drop  become  out  of  phase  with  the 
disturbance,  thereby  tending  to  reduce  the  amplitude  of  the  oscillation.  This  process  repeats  indefinitely  or 
until  the  drop  breaks  up. 

Areas  of  significant  response  occur  for  uj/un  =  near  0.5,  1.0,  2.0,  3.0,  4.0,  and  6.0.  The  peaks  near 
uj/ujn  =  0.5,  1.0,  and  2.0  are  primarily  driven  by  second  mode  oscillations,  while  the  peaks  near  uj/ujn  =3.0, 
4.0,  and  6.0  contain  significant  fourth-mode  effects.  It  is  important  to  note  that  the  fundamental  frequency 
of  the  fourth  mode  is  three  times  that  of  the  fundamental  frequency  of  the  second  mode,  UJ4  =  for  a 
drop  in  a  low  density  gas  (e  <$C  1).  It  is  worth  mentioning  that  a  frequency  ratio  near  5.0  does  not  produce 
a  notable  response. 

For  relatively  low  disturbance  frequencies  [uj  <<  Wn),  a  quasi  steady-state  response  is  observed  whereby 
the  droplet  closely  tracks  the  imposed  perturbation  with  a  very  small  phase  lag.  The  amplitude  (aspect  ratio) 
in  this  region  is  greater  than  the  response  for  higher  off-harmonic  frequency  ratios  due  to  the  negligible  phase 
lag  in  this  region.  At  the  other  extreme  (u;  >>  w^)  the  oscillations  of  the  gas  are  so  fast  that  the  droplet  is 
unable  to  respond  appreciably  to  the  changes  in  the  gas  flow.  Here,  the  droplet  responds  to  a  mean  dynamic 
pressure  generated  by  the  acoustic  wave. 

As  the  frequency  ratio  approaches  a  value  of  0.5,  the  droplet  begins  to  experience  a  sub-harnionic  ex¬ 
citation.  The  frequency  of  its  oscillation  is  0.0604,  which  is  approximately  half  of  the  natural  frequency, 
0.1304.  The  droplet  behavior  is  characterized  by  a  large  deformation  followed  by  a  smaller  one.  The  larger 
response  is  attributed  to  the  peaks  in  the  acoustic  wave,  whereas  the  secondary  response  is  due  to  the  droplet 
attempting  to  maintain  its  natural  frequency.  It  has  been  observed  that  the  second  mode  is  solely  respon¬ 
sible  for  the  activity.  There  is  a  negligible  amount  of  fourth  mode  activity,  arid  the  higher  order  modes  are 
virtually  non-existent  (Murray  1996). 

For  frequency  ratios  near  the  drop  natural  frequency,  the  level  of  droplet  activity  is  high,  such  is  the 
case  for  w/u;^  =  0.79.  Figure  5  displays  the  radial  position  of  the  top  node  for  this  frequency  ratio.  As  the 


64 


droplet  begins  to  deform,  it  becomes  increasingly  excited  as  its  natural  frequency  decreases  and  becomes  in 
tune  with  the  disturbance  frequency.  At  this  point  the  droplet  repsonse  peaks,  and  the  natural  frequency 
continues  to  decrease  and  becomes  out  of  phase  with  the  disturbance.  The  response  begins  to  attenuate, 
which  results  in  the  natural  frequency  becoming  in  tune  with  the  disturbance  again.  The  period  of  this 
envelope  of  oscillations  is  398  seconds  (approximately  6  1/2  oscillations  of  the  imposed  gas  disturbance). 
The  droplet  shapes  at  various  times  in  this  process  are  summarized  in  Figure  6.  Here,  T  is  the  overall  period 
associated  with  the  process  [T  =  398)  and  the  shapes  are  shown  for  various  peaks  and  troughs  in  Fig.  5. 

The  activity  of  the  droplet  is  initially  controlled  by  second  mode  effects.  The  time  history  of  the  mode 
coefficients  for  the  =  0.79  case  is  shown  in  Figure  7.  As  the  oscillatory  amplitude  increaises,  fourth 

mode  coupling  appears.  As  in  all  results  investigated  in  this  study,  6th,  8th,  and  higher  even  modes  gave 
negligible  contributions  to  the  instantaneous  droplet  shapes.  In  conjunction  with  the  nonlinear  frequency 
shift,  the  fourth  mode  effects  have  a  stabilizing  influence  upon  the  droplet.  Figure  8  indicates  that  the 
second  and  fourth  modes  become  coupled  and  destructively  interfere,  which  causes  the  magnitude  of  the 
oscillations  to  dampen  and  the  fourth  mode  effects  to  vanish.  A  similar  envelope  response  is  observed  for  the 
—  0.925  case,  however  the  interference  that  occurs  is  somewhat  constructive,  which  serves  to  amplify 
the  response. 

For  frequency  ratios  between  0.79  and  0.925,  the  second  mode  effects  become  overwhelming  and  the 
droplet  breaks  up.  As  the  frequency  ratio  increases  past  0.925,  the  droplet  activity  drops  off  since  the 
excitation  is  above  the  resonant  condition.  A  region  of  moderate  droplet  response  lies  in  the  frequency  range 
of  1.2  to  1.4,  which  corresponds  to  the  subharmonic  of  the  fourth  mode.  Significant  response  is  also  noted 
at  the  second-harmonic  of  the  fundamental  mode,  as  one  would  expect. 

For  a  frequency  ratio  near  3.0,  the  droplet  experiences  the  third  harmonic  of  the  second  mode  of  oscillation 
as  well  as  the  fundamental  frequency  of  the  fourth  mode  of  oscillation.  The  time  histories  of  the  position  of 
the  top  node  and  the  mode  coefficients  are  presented  in  Figures  9  and  10.  The  overall  period  for  the  process 
is  over  2500  seconds,  during  which  over  150  gas  oscillations  occur.  Throughout  this  process,  the  magnitude 
of  the  second  mode  remains  reasonably  constant,  whereas  the  fourth  mode  grows  and  decays  over  the  period. 
Initially,  the  droplet  oscillates  in  the  second  mode.  As  the  fourth  mode  grows,  it  dominates  the  behavior 
of  the  droplet.  The  fourth  mode  then  decays  to  roughly  the  same  order  of  magnitude  as  the  second  mode. 
Similar  behavior  is  observed  for  a;/u;„  =  3.975. 

Probably  the  most  surprising  result  from  the  frequency  spectrum  is  that  more  response  is  noted  near 
the  second  harmonic  of  the  fourth  mode  ^  6)  than  at  the  primary  fourth  mode  excitation  frequency. 

These  two  peaks  were  the  subject  of  substantial  scrutiny;  model  results  were  replicated  for  several  different 
timesteps  and  mesh  sizes.  The  top  node  position  history  for  this  case  is  shown  in  Figure  11.  Here,  the  overall 
process  takes  about  780  dimensionless  seconds,  corresponding  to  about  95  periods  of  acoustic  excitation.  To 
seek  an  explanation  for  the  heightened  response  (as  compared  to  the  w/wn  =  2.97  case),  the  second  and 
fourth  mode  coefficients  are  compared  in  Figure  12.  As  seen  in  the  upper  plot,  the  second  mode  response  is 
similar  in  both  cases.  However,  the  lower  plot  reveals  a  heightened  fourth  mode  response  near  t  =  400  for  the 
uj/un  =  5.89  case.  Apparently,  there  is  a  constructive  nonlinear  interaction  in  the  uj/ujn  =  5.89  case  leading 
to  stronger  fourth  mode  effects  and  increased  overall  excitation.  Since  viscous  effects  are  more  prevalent  in 
damping  higher  modes,  this  curious  response  may  not  be  replicated  in  actual  experiments. 

4.2  Effect  of  Acoustic  Disturbance  Intensity 

The  influence  of  the  intensity  of  the  acoustic  wave  was  investigated  by  conducting  a  series  of  simulations  at 
fixed  density  and  frequency  ratios,  but  with  varying  Weber  number.  For  these  simulations,  the  time  step 
was  set  to  be  approximately  1/1000  of  the  period  of  a  droplet  oscillation  for  the  m  =  2  mode.  Once  again, 
the  aspect  ratio  was  used  as  an  overall  measure  of  droplet  response;  Figure  13  provides  results  for  the  case 
e  =  0.00123,  to/uJri  =  TO.  During  these  simulations,  it  was  observed  that  droplet  atomization  occurred  at 
Weber  numbers  above  1.10.  Here  it  is  irnportant  to  note  that  this  “critical  Weber  number”  is  dependent 
on  both  density  and  frequency  ratios;  the  1.10  value  is  for  harmonic  excitation  of  a  water  droplet  in  air 
(6  =  0.00123). 

Droplet  break  up  was  studied  by  conducting  trials  using  Weber  numbers  greater  than  the  critical  value. 
From  this  analysis,  three  regimes  of  break  up  were  identified.  Examples  of  these  modes  are  shown  in  Figure 
14.  In  the  “nipple”  breakup  regime  (1.1  <  We  <  2.5),  two  small  satellite  droplets  are  formed  as  a  result  of 
the  nonlinear  motion  of  the  drop.  Since  the  acoustic  intensity  is  barely  above  the  threshold  value,  breakup 
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in  this  region  takes  a  substantial  amount  of  time.  The  nipple  regime  is  similar  to  the  mode  of  break  up 
experienced  in  the  break  up  band  of  frequencies,  as  discussed  in  the  beginning  of  Section  4.1.  As  Weber 
number  is  increased  to  2.5,  atomization  occurs  with  the  formation  of  two  satellite  droplets  which  are  larger 
than  the  central  drop.  This  "kidney'’  regime  occurs  The  kidney  regime  exists  for  2.5  <  We  <  3.0.  Kidney 
regime  breakups  are  generated  in  1-2  oscillations  of  the  acoustic  wave. 

As  the  response  amplifies,  the  drop  flattens  out  in  the  direction  perpendicular  to  the  gas  flow,  and  the 
center  continues  to  flatten  which  results  in  the  pushing  out  of  a  toroid  of  fluid.  Break  up  occurs  when 
the  center  of  the  drop  pinches,  and  the  cross-section  resembles  a  torus.  Atomization  in  this  “toroidal” 
mode  occurs  in  less  than  one  acoustic  period  for  We  >  3.0.  Here,  the  droplet  rapidly  flattens  in  a  plane 
perpendicular  to  the  acoustic  wave.  With  increasing  Weber  number,  the  overall  diameter  of  the  droplet  (at 
the  atomization  point)  increases,  while  the  inner  diameter  of  the  torus  decreases  as  shown  in  Fig.  14.  In 
this  figure,  the  We  =  5.78  case  is  at  a  reduced  scale  for  display  purposes.  The  droplet  shapes  at  high  We 
values  are  consistent  with  “aerodynamic  shattering”  which  has  been  documented  by  observing  the  response 
of  a  droplet  to  a  shock  wave  (see  Hsiang  and  Faeth  (1992)  for  background). 

A  summary  of  the  breakup  times  (nondimensionalized  by  the  period  of  the  acoustic  oscillation,  r)  are 
provided  in  Figure  15.  The  curve  is  somewhat  discontinuous  in  places  due  to  the  fact  that  breakups  occur 
within  discrete  parts  of  a  given  cycle.  Once  again,  for  the  conditions  noted,  no  breakups  occurred  for 
We  <  1.1. 

4.3  Density  Ratio  Effects 

The  effect  of  the  gas/liquid  density  ratio  upon  the  droplet  behavior  was  studied  by  varying  this  parameter 
while  setting  We  =  0.5779  and  —  KO.  At  a  fixed  Weber  number,  the  gas-phase  dynamic  pressure  can 
be  thought  to  be  fixed,  so  that  an  increase  in  density  ratio  under  this  constraint  can  be  viewed  as  a  decrease 
in  density  of  the  droplet.  In  effect,  this  approach  holds  the  dynamic  pressure  constant,  while  decreasing  the 
liquid  inertia.  The  level  of  droplet  response  for  a  range  of  density  ratios  is  shown  in  Figure  16.  As  expected, 
droplet  response  increases  with  density  ratio.  The  results  are  somewhat  discontinuous  as  droplet  response 
enters  the  highly  nonlinear  region,  primarily  because  inflections  in  the  surface  can  start  to  occur  in  this  area. 
For  these  flow  conditions,  droplets  began  to  atomize  for  density  ratios  greater  than  0.1. 

5  Conclusions 

The  model  presented  in  this  paper  has  been  used  to  study  the  nonlinear  evolution  of  an  acoustically  excited 
droplet.  While  the  maximum  response  of  the  droplet  occurs  at  the  harmonic  condition  (u;  =  w^),  significant 
responses  also  occur  for  near  values  of  0.5,  2,  3,  4,  and  6.  The  actual  peak  responses  are  also  slightly 

less  than  these  values  as  a  result  of  the  reduction  in  the  droplet’s  natural  frequency  at  finite  deformation 
amplitude  (nonlinear  frequency  shift).  Droplet  responses  for  the  frequency  ratios  near  0.5  and  2.0  are 
dominated  by  the  second  mode,  whereas  the  coupling  of  the  second  and  fourth  modes  are  responsible  for 
the  responses  to  the  frequency  ratios  near  3.0  and  4.0.  The  droplet  response  for  a;/a;„  «  6,  which  represents 
the  second  harmonic  of  the  fourth  mode,  is  dominated  by  fourth  mode  activity  and  the  amplitude  of  the 
response  actually  exceeds  that  of  some  of  the  lower-order  harmonics. 

As  the  magnitude  of  the  acoustic  disturbance  increases,  droplet  atomization  is  predicted  to  occur.  For 
harmonic  forcing  of  a  liquid  droplet  in  air,  a  critical  Weber  number  of  1.10  divides  the  oscillatory  and 
atomization  regimes.  Droplet  breakups  occurred  in  “nipple”,  “kidney”,  and  “toroidal”  modes  as  Weber 
number  was  increased  above  this  threshold  value.  Break  up  times  were  roughly  inversely  proportional  to 
Weber  number  under  these  conditions.  Finally,  increasing  gas/liquid  density  ratio  under  fixed  Weber  number 
and  frequency  ratio  conditions  was  also  shown  to  heighten  droplet  response. 
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Nomenclature 

a  =■  undisturbed  droplet  radius 
G  —  free  space  Green’s  function 
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h  —  unit  normal  vector 
P  =  pressure 

q  =  velocity  normal  to  local  surface 
t  =  time 

U  =  maximum  velocity  in  acoustic  perturbation 
VVe  =  Weber  number,  \Ne—  pgU’^ a/ cr 
z  =  axial  location 
r  —  radial  location 

a  =  boundary  point  singularity  contribution 

3  =  surface  slope 

r  =  domain  boundary 

€  =  ratio  of  gas  to  liquid  density 

K  =  surface  curvature 

p  =  density 

(T  =  surface  tension 

<j)  ~  velocity  potential 

ujri  =  undisturbed  droplet  second  mode  natural  frequency 
ijj  =  acoustic  frequency  (gas  phase) 

Subscripts 

g  =  gas  phase  properties 
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1.  Schematic  of  Computational  Domain  and  Boundary  Conditions. 

2.  Comparisons  of  Predicted  Frequency  Shifts  for  Nonlinear  Oscillations  of  a  Droplet. 

3.  Comparisons  of  Steady  State  Droplet  Profiles. 

4.  Frequency  Response  Spectrum. 

5.  Time  History  of  Top  Node  Position. 

6.  Time  History  of  Droplet  Profiles. 

7.  Time  History  of  Mode  Coefficients. 

8.  Oscillatory  Mode  Interference. 

9.  Time  History  of  Top  Node  Position. 

10.  Time  History  of  Mode  Coefficients. 

11.  Time  History  of  Top  Node  Position. 

12.  Second  and  Fourth  Mode  Coefficients. 

13.  Effects  of  the  Weber  Number  upon  Droplet  Behavior  (e  =  0.00123,  w  =  w^i). 

14.  Weber  Number  Induced  Droplet  Break  Up  Profiles  (e  =  0.00123, w  =  u;n). 

15.  Time  Required  for  Droplet  Break  Up  for  Various  Weber  Numbers  (e  =  0.00123, u; 

16.  Influence  of  Liquid  Density  upon  Droplet  Behavior.  (We  =  0.5779,  cj  =  u;^). 
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Figure  1:  Schematic  of  Computational  Domain  and  Boundary  Conditions 


Figure  2:  Comparisons  of  Predicted  Frequency  Shifts  for  Nonlinear  Oscillations  of  a  Droplet 
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Figure  9:  Time  History  of  Top  Node  Position 
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Figure  10:  Time  History  of  Mode  Coefficients 
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Position  of  Top  Node 


^  =  5.89,  We  =  0,5779,  £  =  0.00123 


Figure  11:  Time  History  of  Top  Node  Position 
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Figure  12:  Second  and  Fourth  Mode  Coefficients 
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10  Appendix  E  -  BEMs  for  Two-Fluid  Flows 


Heister,  S.  D.,  “Boundary  Element  Methods  for  Two-Fluid  Free  Surface  Flows”,  In  Re¬ 
view,  Engineering  Analysis  with  Boundary  Elements^  1996. 
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Abstract 

Flows  in  which  a  low  density  (gaseous)  phase  contributes  to  the  motion  of  a  free  surface  have  been  studied 
numerically  using  boundary  element  methods.  A  stable,  time-accurate  integration  scheme  for  the  coupled, 
nonlinear  free  surface  boundary  conditions  is  described  for  the  case  where  both  phases  are  incompressible  and 
inviscid.  The  effect  of  increasing  gas/liquid  density  ratio  on  stability  of  the  methodology  is  briefly  discussed. 
The  technique  is  illustrated  through  a  series  of  examples  involving:  an  infinite-length  (periodic)  liquid  jet 
in  a  “wind-induced”  flow  regime;  a  two-dimensional  liquid  column  subjected  to  acoustic  excitation;  and  a 
finite-length  liquid  jet  injected  into  a  quiescent  gas. 


Introduction 

Within  the  past  two  decades,  researchers  have  adopted  boundary  element  solutions  to  problems  involving 
nonlinear  deformations  of  a  free  surface.  In  these  problems,  the  boundary  element  approach  is  attractive 
because  the  “grid”  is  simplified  by  one  dimension  over  that  of  a  more  traditional  computational  fluid  dynamic 
calculation.  For  two-dimensional  geometries,  the  grid  in  the  BEM  formulation  is  simply  a  curved  line.  In 
problems  where  surface  topology  changes  occur  (such  as  atomization),  this  simplification  becomes  quite 
valuable.  In  fact,  recent  BEM  simulations^  have  been  conducted  in  which  calculations  proceed  beyond 
atomization  events. 

Numerical  approaches  utilizing  Finite  Element  Methods  (FEM)  and  the  Volume  of  Fluid  (VOF)  method 
have  also  been  applied  to  problems  of  this  type.  The  VOF  technique  relies  on  the  interpolation  of  the  surface 
location  from  a  fixed  computational  mesh.  In  capillary  flows,  this  interpolation  procedure  can  introduce 
substantial  inaccuracy  in  determining  surface  curvature  (and  hence  the  capillary  force).  As  an  example,  a 
typical  VOF  calculation  involving  sloshing  of  a  fluid  in  a  tank“  exhibits  a  1%  error  in  preserving  the  liquid 
volume,  while  a  BEM  calculation  exhibits  a  0.01%  volume  error^.  A  typical  BEM  calculation^  involving 
nonlinear  oscillations  of  a  droplet  utilizing  45  nodes  achieves  a  maximum  volume  error  of  0.04  %,  while  a 
comparable  FEM  calculation^  would  use  over  1600  nodes  and  produce  a  volume  error  of  0.8%. 

For  these  reasons,  we  have  seen  a  variety  of  applications  of  BEMs  to  problems  involving  large  deformations 
of  a  free  surface.  Several  solutions  have  been  developed  for  nonlinear  evolution  of  water  waves^^^,  and  for 
nonlinear  deformations  of  both  viscous  and  inviscid  drops^’®.  A  variety  of  solutions  have  been  obtained  for 
creeping  (Stokes)  flows  in  liquid  columns^" and  in  annular  layers^ Inviscid  solutions  have  also  been 
obtained  for  both  infinite^ ^  and  finite-length^  liquid  jets,  as  well  as  for  dripping  flows^  fountains^  and  fluid 
sloshing  problems^^. 

These  simulations  have  all  had  to  treat  the  primary  nonlinearity  associated  with  free  surface  flows  involv¬ 
ing  the  fact  that  the  current  surface  location  cannot  be  decoupled  from  the  pressure/ velocity  field.  In  other 
words,  the  surface  pressure  is  dependent  on  the  shape  of  the  interface,  which  in  turn  depends  on  the  local 
pressure  distribution.  Most  researchers  have  treated  this  nonlinearity  by  utilizing  small  time  steps  such  that 
the  surface  pressure/ velocity  is  invariant  over  a  given  step.  More  elaborate  treatments  have  been  suggested 
and  validated  by  Liggett  and  coworkers^^’^^. 

Recently,  a  variety  of  models  have  been  developed  in  order  to  include  the  influence  of  gas-phase  pressure 
distribution  on  the  nonlinear  evolution  of  the  interface^^”^®.  Developing  a  capability  to  address  these  flows 
permits  the  consideration  of  problems  in  which  wind-induced,  or  acoustic  interactions  from  the  gas  phase 
can  be  included  as  physical  factors  affecting  distortion  of  the  interface.  In  this  case,  the  nonlinear  free 
surface  boundary  condition  is  coupled  between  the  two  fluids  bordering  the  interface.  This  complication 
introduces  some  unique  computational  issues  for  those  interested  in  modeling  these  flows.  In  this  paper,  we 
address  these  issues  as  applied  to  several  problems  of  engineering  interest.  The  following  section  provides 
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a  detailed  description  of  the  BEM,  while  subsequent  sections  highlight  boundary  condition  treatments  for 
three  example  problems. 


Model  Development 


Our  interest  lies  in  developing  models  which  can  address  capillary  (surface  tension)  forces  at  the  interface. 
By  choosing  the  liquid  density  (p),  the  size/height  of  the  liquid  body  (a),  and  a  farfield  velocity  {[/)  as 
dimensions,  the  gas/liquid  density  ratio: 

^  =  Pg/P  (1) 

and  the  Weber  number  based  on  gas  density: 


(2) 


become  the  two  dimensionless  parameters  characterizing  the  flowfield.  Here,  the  Weber  number  measures  the 
ratio  of  inertial  forces  imposed  by  the  gas  phase  to  surface  tension  ((j)  forces.  In  the  following  development, 
we  presume  that  the  nondimensionalization  described  above  has  been  applied. 

We  assume  that  both  phases  can  be  represented  as  incompressible,  inviscid  fluids.  In  this  case,  a  velocity 
potential  (whose  gradient  is  simply  the  velocity)  exists.  Let  <t>  and  (f>g  represent  the  velocity  potential  in 
liquid  and  gaseous  phases,  respectively.  Continuity  requires  that  both  velocity  potentials  satisfy  Laplace  s 
equation: 

V-0  =:  =  0  (3) 

Proceeding  with  a  standard  BEM  formulation,  the  integral  form  of  eqn  3  for  the  liquid  domain  becomes: 


^-,G]dr  =  o 


(4) 


where  0(n)  is  the  value  of  the  potential  at  a  point  T  denotes  the  boundary  of  the  domain,  and  G  is  the 
free-space  Green’s  function  corresponding  to  Laplace’s  equation.  An  analogous  form  of  eqn  4  can  also  be 
derived  for  the  gas  phase  potential.  For  a  well-posed  problem,  either  <i>  oi  q  —  d<t>ldn  must  be  specified  at 
each  “node”  on  the  boundary.  Here  n  is  the  outward  normal  to  the  boundary  so  that  q  represents  the  velocity 
normal  to  the  boundary.  The  quantity  a  in  eqn  4  results  from  singularities  introduced  as  the  integration 
passes  over  the  boundary  point,  fl. 

Models  have  been  developed  for  both  two-dimensional  and  axisymmetric  flowfields.  If  we  let  r  and  z 
denote  radial  and  axial  coordinates,  respectively,  and  denote  the  base  point  with  subscript  “i”,  the  Green’s 
function  solution  to  the  axisymmetric  Laplacian  can  be  written: 


ArK{p) 

y/{r  +  +  {z  -  Zi)^ 


where 


(r- +  r,)2  +  (z  -  z,)2 


(6) 


and  K{p)  is  the  complete  elliptic  integral  of  the  first  kind.  For  computational  efficiency,  this  quantity  is 
calculated  using  a  curve  fit^^  which  has  an  accuracy  to  10"’®. 

In  the  case  of  a  2-D  flow  (letting  x  and  y  represent  the  coordinates),  we  have: 


G  =  ’■^1  =  ~  4-  {y 


(7) 


In  both  c2Lses,  we  have  utilized  linear  elements  in  the  formulation  of  a  set  of  integral  equations  over 
a  discretized  boundary.  Linear  elements  are  desirable  for  problems  with  moving  surfaces  since  the  BEM 
solution  in  this  case  will  return  velocities  at  the  ends  of  each  segment  corresponding  to  the  nodal  locations. 
In  the  case  of  constant  elements  (for  example),  a  nodal  velocity  would  be  calculated  at  the  center  of  each 
segment,  but  the  end  points  of  the  segment  would  require  a  separate  treatment  (based  on  interpolation  of 
nodal  velocities)  to  update  the  mesh  at  subsequent  times. 
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For  the  2-D  flows,  integration  across  a  segment  can  be  carried  out  analytically.  Singularities  resulting  from 
integration  across  a  segment  containing  the  base  point  are  also  integrable.  In  the  case  of  axisymmetric  flow, 
the  integrations  must  be  carried  out  numerically.  In  this  case,  we  choose  a  four-point  Gaussian  quadrature 
for  evaluation  of  integrals.  Logarithmic  singularities  which  arise  in  the  elliptic  integral  when  the  segment 
contains  the  base  point  are  treated  with  a  special  Gaussian  integration  designed  to  accurately  treat  this 
condition.  Additional  details  regarding  the  numerical  implementation  can  be  found  in  Refs,  1  and  17  for 
the  axisymmetric  and  2-D  cases,  respectively. 


Free  Surface  Treatment 


Governing  Equations 

The  main  challenge  in  developing  models  capable  of  tracking  large  deformations  of  an  interface  lies  in  the 
treatment  of  the  free  surface  itself.  Since  capillary  forces  are  important,  it  is  crucial  to  develop  a  treatment 
capable  of  accurately  determining  surface  curvature  at  all  times  during  the  simulation.  For  this  reason,  all 
models  employ  fourth-order  centered  differencing  (on  a  generalized,  variable  spacing  mesh)  to  determine 
surface  curvature.  Curvature  is  calculated  based  on  coordinate  derivatives  as  a  function  of  distance  along 
the  surface  using  the  parametric  representation  due  to  Smirnov“®. 

The  modeler  has  the  choice  of  tracking  the  motion  of  free  surface  nodes  in  a  variety  of  directions®.  In 
current  models,  we  have  opted  to  track  surface  nodes  along  the  local  liquid  velocity  vector.  Under  this 
assumption,  for  an  axisymmetric  situation,  flow  kinematics  require: 

Dt  dz  Dt  dr 

where  the  notation  D{)/ Di  implies  a  Lagrangian  derivative  for  points  on  the  surface  moving  with  the  local 
liquid  velocity. 

Recognizing  that  our  BEM  solver  will  return  velocities  normal  to  the  surface,  we  employ  the  velocity 
transformations: 

^  “  qsin{l3)  (9) 


where  (3  is  the  local  wave  slope  and  d<l)/ds  is  the  velocity  tangential  to  the  local  surface.  This  tangential 
velocity  is  calculated  using  5-point  centered  differences  on  <f>,  except  for  nodes  adjacent  to  ends  of  the  free 
surface,  where  a  3-point  formula  is  employed.  The  local  wave  slope,  (3,  is  calculated  following  the  formulation 
of  Medina^^.  For  each  node,  a  parabola  is  defined  such  that  it  passes  through  the  previous  node,  the  node 
in  question,  and  the  following  node.  The  slope  of  the  surface  is  given  by  the  tangent  to  the  parabola  at  the 
central  node.  A  completely  analagous  treatment  can  be  employed  for  2-D  flows  by  replacing  r  with  y  and  2 
with  X  in  eqns  8  and  9. 

Dynamics  of  the  interface  are  addressed  through  the  unsteady  Bernoulli  equation.  In  an  Eulerian  system 
where  time  derivatives  are  assumed  to  occur  at  a  fixed  spatial  location,  the  dimensionless  form  of  this  relation 
for  the  liquid  surface  is: 


(10) 


where  Pg  is  the  ga^  pressure  at  the  interface,  and  k  is  the  surface  curvature.  The  Eulerian  -  Lagrangian 
transformation  for  nodes  on  the  interface  moving  with  the  liquid  velocity  is: 


Dt  dt 


+  •  V(  ) 


Using  this  transformation,  the  Bernoulli  equation  in  the  liquid  becomes: 


D<f) 

~Di 


while  an  analagous  treatment  for  the  gas  phase  gives: 

-  P, 


(11) 


(12) 


(13) 
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Time  Integration  Scheme 

Mathematically,  eqns  8,  12,  and  13  provide  a  system  of  relations  to  describe  the  evolution  of  the  surface  shape 
(r,  ^  or  y,z)  and  velocity  potentials  for  unsteady  motion  of  the  interface.  These  equations  are  integrated 
in  time  using  a  fourth-order  Runge-Kutta  schemed  This  scheme  has  the  advantage  of  full  fourth-order 
accuracy  without  the  requirement  of  a  knowledge  of  the  “history”  for  a  given  nodal  location;  i.e.  information 
at  previous  time  levels  is  not  required  in  the  integration  algorithm.  This  feature  can  be  advantageous  in 
calculations  where  a  variable  timestep  is  employed  or  where  the  number  of  nodes  along  the  free  surface  is 
not  constant  (due  to  atomization  events  or  surface  regridding). 

As  mentioned  previously,  the  main  challenge  in  this  problem  is  the  development  of  a  stable,  consistent 
procedure  to  handle  the  coupled,  nonlinear  boundary  conditions  at  the  interface  (eqns  12  and  13).  More 
specifically,  if  we  regard  eqn  13  as  an  expression  for  Pg,  than  an  approximation  for  the  derivative  D(f>g/Dt 
is  required.  We  have  found,  that  for  a  wide  array  of  problems  it  is  adequate  to  approximate  this  derivative 
using  a  first-order  backward  difference  scheme. 


D<p,  _ 

Dt  At 


(14) 


where  ‘T’  denotes  time  level. 

In  problems  where  the  gas  flows  tangential  to  the  liquid  surface  over  most  of  the  time-dependent  process, 
the  major  term  influencing  the  gas  pressure  in  eqn  13  turns  out  to  be  the  term  involving  {V(f>g)-,  In  this  case, 
highly  accurate  approximation  of  D<f)gf  Dt  is  obviously  not  required.  However,  we  have  found  the  first-order 
backward  difference  (eqn  14)  to  be  adequate  even  in  problems  where  stagnation  points  are  present.  In  this 
case,  inspection  of  eqn  13  indicates  that  the  D(t)g/Dt  term  is  the  major  component  leading  to  changes  in  gas 
pressure.  In  the  examples  which  follow,  we  show  that  the  simple  approximation  in  eqn  14  is  adequate  for  a 
wide  variety  of  flows. 

The  following  procedure  is  implemented  for  nodes  on  the  free  surface: 

•  At  the  start  of  a  given  time  step,  the  value  of  (j>  is  known.  Using  this  value  as  a  boundary  condition 
on  the  interface,  the  liquid  velocity  q  can  be  determined  via  solution  of  Laplace’s  equation  (eqn  4). 

♦  Since  the  gas  nodes  on  the  interface  are  fixed  to  move  with  the  liquid  nodes,  this  liquid  velocity  is 
used  as  the  gas  phase  boundary  condition  {qg  =  -q)  to  calculate  the  (j)g  value  on  the  gas  side  of  the 
interface. 


•  This  value  of  the  gas  phase  velocity  potential  is  then  used  in  eqn  13  to  determine  the  gas  pressure  at 
this  new  time  step  using  the  approximation  for  D4>glDi  given  in  eqn  14. 

•  The  gas  pressure  at  the  new  time  is  then  used  in  eqn  12  to  calculate  the  current  D^jDi  which  is  then 
integrated  in  time. 


Since  the  nodes  on  the  interface  are  allowed  to  move  with  their  local  velocity,  over  time  they  tend  to  group 
themselves  in  regions  of  high  curvature.  This  phenomena  leaves  regions  of  lower  curvature  poorly  defined.  To 
alleviate  this  problem,  the  surface  mesh  is  regridded  using  a  series  of  cubic  splines  (for  surface  coordinates, 
0,  and  <i>g)  at  each  time  step  to  keep  the  spacing  between  the  nodes  constant  along  the  surface.  The  use 
of  the  Runge-Kutta  integration  scheme  is  well  suited  to  this  type  of  remeshing,  since  it  does  not  require 
information  on  node  positions  at  previous  time  levels  to  predict  the  subsequent  motion  of  the  surface.  Also, 
we  note  that  regridding  the  surface  can  be  accomplished  in  this  case  since  the  approximation  for  D(j)g/  Dt 
(eqn  14)  involves  only  two  time  levels.  If  more  accurate  representations  of  this  derivative  are  required,  then 
regridding  tends  to  destroy  information  about  previous  <j>g  values  on  given  nodal  locations. 

Finally,  we  note  that  the  regridding  process  does  provide  a  natural  ^‘smoothing”  of  the  surface.  Many 
previous  authors^’®’^^*^^  have  been  forced  to  implement  smoothing  procedures  to  alleviate  “zig-zag”  instabili¬ 
ties  which  develop  on  the  surface  after  a  large  number  of  time  steps.  Because  of  the  regridding  procedure,  we 
have  not  had  to  implement  any  formal  smoothing  of  the  surface  (or  any  other  functions  associated  with  the 
surface).  For  these  reasons,  calculations  using  the  methodology  described  above  have  very  little  numerical 
dissipation.  In  the  following  sections,  we  provide  three  examples  to  illustrate  the  results  of  the  free  surface 
treatment  described  above. 
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Example  #1:  Infinite  Liquid  Jet  in  Wind-Induced  Flow  Regime 

The  distortion  and  atomization  of  a  liquid  jet  is  one  of  the  most  frequently  studied  free  surface  problems 
due  to  its  wide  array  of  applications  (atomizers,  ink-jet  printers,  fire  hoses,  etc.).  To  date,  most  analyses 
have  presumed  a  periodic  solution  exists,  whereby  the  study  focuses  on  a  single  wavelength  of  fluid.  Such 
conditions  can  be  experimentally  simulated  in  the  case  of  low  speed  jets  by  providing  low  amplitude  acoustic 
or  mechanical  excitation  of  the  orifice““.  Until  recently,  the  nonlinear  analyses  applied  to  this  problem 
have  ignored  the  presence  of  the  gas  phase.  Using  the  formulation  described  in  the  previous  sections,  we 
have  developed  a  fully  nonlinear  simulation  of  the  jet  in  the  presence  of  a  gas  (typically  referred  to  as  the 
wind-induced  flow  regime). 

A  schematic  of  the  computational  domain  is  shown  in  Fig.  1  in  which  liquid  nodes  are  denoted  “o”,  while 
the  gas  nodes  are  denoted  ‘V’.  We  presume  that  the  coordinate  system  is  attached  to  the  liquid,  so  that 
U  is  the  relative  velocity  between  the  two  fluids  at  a  location  far  from  the  interface.  We  place  the  upper 
boundary  in  the  gas  phase  far  enough  away  from  the  surface  such  that  there  is  negligible  flux  of  fluid  through 
the  interface.  The  assumed  periodicity  requires  that  on  the  vertical  boundaries  gr  =  0  in  the  liquid  phase  and 
ddgidr  =:  0  in  the  gas  phase.  No  nodes  are  required  along  the  centerline  since  both  G  and  dGjdn  vanish 
for  the  axisymmetric  geometry  assumed. 

The  initial  wave  shape  of  the  surface  is  assumed  to  be: 


r  =  1  -f  T]cos{kz) 


(15) 


where  t]  is  the  initial  surface  deflection  and  k  is  the  wave  number  of  the  forcing  disturbance  on  the  jet.  For 
7/  <<  1,  results  from  a  linear  analysis^^  can  be  employed  to  give  the  initial  gas-phase  pressure  and  velocity 
potential  along  the  interface: 
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where  A'o  and  Ki  are  modified  Bessel  functions  of  the  second  kind.  In  the  liquid  phase,  we  presume  that 
the  interface  is  initially  at  rest  to  provide  a  starting  condition  for  the  integration. 

Typical  calculations  involve  60-100  nodes  along  the  interface.  Vertical  surfaces  in  the  gas  phase  employ 
exponential  stretching  to  permit  greater  resolution  near  the  interface.  Near  the  interface,  nodal  spacings 
along  vertical  surfaces  are  kept  at  values  near  that  of  the  interface  itself.  A  dynamic  time-step  criteria  is 
employed  to  insure  that  nodes  move  no  more  than  5%  of  the  nodal  spacing  on  a  given  time  step.  This 
criteria  permits  considerable  acceleration  of  the  code  (as  compared  to  the  constant  time  step  case)  since  the 
initial  surface  velocities  are  very  small.  Surface  velocities  grow  to  very  large  values  as  curvature  increases 
dramatically  near  the  pinching  of  the  jet.  Calculations  are  stopped  when  a  node  reaches  a  distance  of  1%  of 
the  initial  jet  radius  from  the  centerline.  A  typical  calculation  takes  about  25,000  CPU  seconds  on  an  IBM 
RISC  6000  Model  580  machine. 

Results  from  sample  calculations  are  included  in  Figs  2  and  3.  In  Fig.  2,  the  evolution  of  the  jet  is  shown 
for  We  =  1,  e  =  0.00129  which  would  correspond  to  a  1  mm  water  jet  issuing  at  7.7  m/s  into  ambient  pressure 
air  (multiple  waves  are  shown  here  for  clarity).  A  dimensionless  wave  number  of  k  =  1.027  is  selected,  since 
linear  theory^^  predicts  this  to  be  the  most  unstable  value  for  this  flow  condition.  This  calculation  employs 
65  nodes  along  the  interface  and  a  total  of  87  and  259  nodes  in  liquid  and  gaseous  domains,  respectively. 
Under  these  conditions,  aerodynamic  and  surface  tension  forces  are  comparable  and  the  simulation  predicts 
that  the  jet  will  break  at  the  centerline  leading  to  the  formation  of  a  “main”  and  “satellite”  drop  from 
each  wave  along  the  surface.  Jets  atomizing  in  this  fashion  are  said  to  lie  in  the  first  wind-induced  regime. 
Predictions  of  the  size  of  main  and  satellite  drops  agree  well  with  experimental  measurements^®  for  low  speed 
jets. 

Figure  3  shows  the  influence  of  increased  gas  velocity  for  the  case:  We  =  5,  c  =  0.0013,  k  =  4.23.  In 
this  case,  the  increased  influence  of  the  gas-phase  leads  to  the  pinching  of  an  annular  ring  of  fluid  at  the 
periphery  of  the  jet  (characteristic  of  the  second  wind-induced  regime).  The  increased  complexity  of  the 
surface  in  this  case  necessitated  a  more  refined  grid  which  employs  89  nodes  along  the  interface.  While 
viscous  effects  would  tend  to  cause  the  annular  ring  to  be  swept  downstream,  aerodynamic  ripples  of  this 
nature  have  been  observed  experimentally^®. 


84 


Spangler^^  has  investigated  the  influence  of  e  on  the  performance  of  the  code.  The  algorithm  produces 
stable  results  for  c  values  approaching  unity,  but  diverges  rapidly  for  values  near  or  above  this  threshold.  As 
e  increases,  the  importance  of  the  gas  phase  also  increases.  For  e  >  1  we  presume  it  would  be  more  logical  to 
interchange  the  treatment  of  liquid  and  gas  phases  in  the  procedure  described  above.  In  addition,  it  is  likely 
that  more  accurate  and  stable  results  would  be  obtained  by  tracking  nodes  along  the  '‘gas”  phase  velocity 
vector  for  these  conditions.  In  practice,  e  values  less  than  0.01  characterize  nearly  all  gas/liquid  flows;  even 
at  high  pressure  conditions  present  in  a  rocket  combustion  chamber. 


Example  ^2:  Liquid  Column  Subjected  to  Acoustic  Perturba¬ 
tions 

Nonlinear  interactions  of  various  atomization  processes  with  acoustic  waves  have  been  implicated  as  a  possible 
source  of  combustion  instabilities  in  liquid  rocket  engines-^.  For  this  reason,  a  2-D  model  was  developed 
to  assess  the  distortion  of  a  liquid  column  subjected  to  a  transverse  acoustic  wave.  Here,  we  presume  that 
the  wavelength  of  the  acoustic  disturbance  is  much  greater  than  the  diameter  of  the  column  such  that  the 
acoustic  wave  can  be  modeled  as  an  unsteady  incompressible  flow. 

Under  these  assumptions,  the  computational  domain  for  this  analysis  is  shown  in  Fig.  4.  Due  to  symme¬ 
try,  only  a  1/4  cylinder  domain  need  be  considered.  In  Fig.  4,  liquid  nodes  are  represented  by  “o”  and  gas 
nodes  by  For  this  2-D  problem  (flow  is  into  the  page),  nodes  are  required  along  both  symmetry  axes. 
The  boundary  condition  ^  =  0  is  applied  along  the  y  =  0  symmetry  plane.  Along  the  a:  =  0  plane,  we 

specify  ^  =  0  and  d(f)gldy  =  0  for  the  gas  phase.  Once  again,  exponential  stretching  is  applied  to  gas  nodes 
along  the  symmetry  axis  in  order  to  increase  resolution  near  the  stagnation  point  at  the  junction  with  the 
column. 

The  outer  boundary  in  the  gas  domain  is  placed  far  enough  (15  jet  radii)  from  the  liquid  such  that  the 
velocity  potential  on  this  surface  may  be  assumed  to  be  a  pure  acoustic  wave  traveling  in  the  horizontal 
direction: 

(i)g  =  xcos{ujgt/2)  (18) 

where  ujg  is  the  frequency  of  the  acoustic  wave.  The  factor  of  1/2  in  the  cosine  wave  results  from  the  fact 
that  the  solution  is  insensitive  to  flow  direction,  i.e.  a  wave  traveling  to  the  right  sets  up  the  same  pressure 
distribution  as  a  wave  traveling  to  the  left.  In  this  problem,  we  expect  large  column  response  when  the 
frequency  of  the  acoustic  wave  is  near  the  natural  frequency  of  the  column.  A  linear  analysis^^  shows  that 
the  dimensionless  column  natural  frequency  can  be  expressed: 


I 

\l{l-\-€)We 


(19) 


Here  we  have  chosen  the  column  radius,  liquid  density,  and  the  peak  velocity  from  the  acoustic  wave  as 
dimensions. 

A  time  step  of  0.005  was  used  in  all  simulations.  The  model  was  validated  using  the  analytic  solution  for 
flow  over  a  cylinder.  In  addition,  results  compare  very  well  with  the  steady-state  solutions^^  for  the  case  of 
a  uniform  flow  over  a  liquid  cylinder.  Typical  grids  employ  25  nodes  on  the  interface,  21  nodes  on  the  outer 
gas  boundary,  15  gas  nodes  along  the  symmetry  planes,  and  8  liquid  nodes  along  symmetry  planes.  Typical 
run  times  are  4-16  hours  on  a  Sun  Sparcstation  1000  computer. 

Sample  results  from  the  model  are  provided  in  Figs.  5-7.  In  Fig.  5,  the  motion  of  a  node  at  the  top  of 
the  column  is  shown  for  the  case:  We  —  0.1,  e  =  0.01,  Ug  =  u;„  which  represents  excitation  of  the  column 
at  its  lowest-order  natural  frequency.  Even  though  the  energy  input  is  unbounded,  the  column  response 
is  bounded.  This  phenomena  is  caused  by  a  shift  in  the  column  natural  frequency  at  larger  deformations. 
At  large  amplitudes,  the  natural  frequency  is  substantially  lower  than  that  predicted  by  the  linear  theory 
(eqn  19),  which  ultimately  leads  to  a  destructive  interference  from  the  constant  frequency  acoustic  wave. 
Nonlinear  frequency  shifts  have  been  characterized  for  both  droplets^  and  columns^^  in  the  literature.  The 
simulation  fails  at  ^  «  95  due  to  the  development  of  a  sawtooth  appearance  along  the  interface.  Such 
occurances  are  common  in  long  time  integrations  of  this  nature  (note  that  we  have  taken  almost  20,000  time 
steps  prior  to  failure  of  the  simulation).  As  mentioned  previously,  numerical  smoothing  would  be  required 
to  extend  the  simulation  to  longer  times. 
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The  shapes  of  the  column  at  various  stages  in  the  process  are  shown  in  Fig.  6.  In  contrast  to  oscillations  of 
droplets,  the  column  is  fairly  insensitive  to  fourth-mode  activity  so  that  a  second  (primary)  mode  oscillation 
occurs  at  all  gas  frequencies^'.  The  effects  of  gas  density  changes  are  assessed  in  Fig.  7  for  the  harmonic 
case  (ijg  =  by  comparing  the  maximum  aspect  ratio  of  the  column  {AR{max))  at  various  e  values.  Since 
Weber  number  is  fixed  in  these  simulations,  the  effects  of  increasing  e  can  be  thought  of  as  a  decreased  liquid 
density  (and  hence  liquid  inertia).  For  this  reason,  we  expect  a  greater  column  response  at  higher  e  values. 
Figure  7  confirms  this  behavior,  showing  a  greater  sensitivity  to  €  as  c  —  0.1.  At  c  >  0.01,  the  treatment  of 
the  interface  does  not  provide  viable  solutions  and  numerical  zig-zag  instabilities  are  generated  early  in  the 
calculation. 

It  is  interesting  to  note  that  these  solutions  break  down  at  a  much  lower  value  of  e  than  those  in  the  first 
example.  While  the  breakdown  point  is  clearly  dependent  on  Weber  number,  another  explanation  for  this 
behavior  is  the  increased  sensitivity  of  the  solution  to  the  approximation  of  D(j)g/Dt  near  the  stagnation 
point  at  the  base  of  the  column.  In  fact,  Rutz^'  has  extended  We  =  0.1  solutions  to  higher  e  values  by 
removing  the  regrid  procedure  and  implementing  fourth-order  time  integrations  of  eq.  13. 


Example  #3:  Liquid  Jet  Injected  into  Quiescent  Gas 


A  final  example  providing  unique  boundary  conditions  for  a  fully-coupled  gas/liquid  flow  is  outlined  in  Fig. 8. 
In  this  axisymmetric  problem,  a  finite-length  liquid  jet  is  injected  into  a  quiescent  gas.  Note  that  the  orifice 
geometry  can  also  be  included  in  this  simulation;  a  unique  treatment  as  compared  to  most  atomization 
models.  The  outer  boundary  in  the  gas  domain  is  made  large  enough  to  encompass  the  entire  liquid  domain 
for  the  entire  duration  of  the  calculation.  In  addition,  this  boundary  is  placed  far  enough  from  the  liquid 
so  it  can  be  assumed  to  be  at  constant  pressure.  Similarly,  the  hemispherical  inflow  boundary  is  placed  far 
enough  from  the  orifice  entrance  such  that  constant  pressure  can  be  assumed. 

Combustion  instabilities  have  been  attributed  to  the  “dynamic  orifice  flow”  created  by  a  time-varying 
discharge  pressure.  Using  the  domain  in  Fig.  8,  we  can  assess  liquid  behavior  for  both  constant  and  time- 
varying  orifice  pressure  drops.  Even  though  a  stagnation  point  lies  at  the  tip  of  the  jet,  we  have  found  that 
the  free-surface  treatment  described  in  Example  #  1  is  adequate  for  this  problem.  Nodes  are  added  along 
the  free  surface  as  the  jet  issues  from  the  orifice  so  as  to  maintain  a  roughly  constant  nodal  spacing.  Surface 
nodes  are  regrid  at  each  time  step  by  fitting  surface  coordinates  and  velocity  potentials  with  cubic  splines, 
as  described  in  Example  #  1. 

The  inflow  and  outer  gas  boundaries  require  a  unique  treatment  in  this  example.  On  the  inflow  boundary, 
Bernoulli’s  equation  can  be  written: 

I  =  i(V«=  -  (20) 

where  Pio  is  the  prescribed  inflow  pressure.  On  the  outer  gas  boundary,  Bernoulli’s  equation  is: 

^  (21) 


where  Pgo  is  the  prescribed  gas  pressure.  Equations  20  and  21  require  no  Eulerian/Lagrangian  transformation 
since  nodes  on  these  boundaries  remain  fixed. 

In  this  example,  we  presume  that  all  pressures  are  nondimensionalized  using  the  liquid  density  and  the 
ideal  orifice  exit  velocity  (ctssuming  a  discharge  coefficient  of  unity).  Under  this  assumption,  Pio  and  Pgo 
must  be  related  by: 


Plo  -  Pgo  =  - 


(22) 


where  RinUt  is  the  radius  of  the  hemispherical  inlet  at  the  entry  to  the  orifice.  If  an  oscillating  Pgo  is  specified, 
then  the  mean  value  is  used  in  eqn  22  to  find  the  appropriate  Pio  to  support  the  desired  nondimensionalization. 
Using  this  approach,  the  problem  is  characterized  by  a  Weber  number,  density  ratio,  and  the  frequency  and 
amplitude  of  any  pressure  perturbation  applied  at  the  outer  gas  boundary. 

On  the  liquid  inflow  boundary,  we  begin  calculations  by  setting  q  =  —1/(272?^^^^)  along  the  inflow 
boundary  and  2issume  a  small  column  of  liquid  (with  a  hemispherical  cap)  is  issuing  from  the  orifice  with 
constant  axial  velocity,  <f>  =  z.  We  set  g  =  0  along  solid  boundaries  representing  walls  of  the  orifice  or  of  the 
chamber  and  solve  Laplace’s  equation  to  provide  initial  values  for  (f)  along  the  inflow  boundary.  The  outer 
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gas  boundary  is  placed  far  enough  from  the  small  liquid  jet  such  that  a  stagnant  condition  (0^  =  0)  initially. 
With  this  information,  we  can  begin  the  time-stepping  procedure  described  in  Example  #1  to  determine  (pg 
values  and  velocities  along  the  interface.  In  this  case,  we  integrate  an  additional  two  equations  (eqn  20  and 
21)  to  update  (p  on  the  inflow  boundary  and  <pg  on  the  outer  gas  boundary.  These  expressions  are  integrated 
using  the  fourth-order  Runge  Kutta  integration  described  previously. 

Typical  calculations  employ  a  grid  spacing  approximately  equal  to  20%  of  the  orifice  radius,  with  a  time 
step  of  0.005.  Run  times  can  be  quite  substantial  (1-2  weeks)  due  to  the  fact  that  the  size  of  the  gas/liquid 
interface  increases  with  time.  It  is  not  uncommon  to  have  grids  with  several  hundred  nodes  in  both  liquid 
and  gas  phases. 

Results  of  sample  calculations  using  the  model  are  presented  in  Figs.  9  and  10.  Figure  9  summarizes 
a  series  of  runs  aimed  at  investigating  the  initial  behavior  of  the  jet  for  a  variety  of  density  ratios.  As  e 
is  increased,  the  momentum  required  to  displace  the  gas  also  increases,  thereby  yielding  the  “mushroom'’ 
shaped  jet  tips  shown  in  Fig.  9.  In  Fig.  10,  the  effect  of  unsteadiness  in  gas  pressure  is  investigated.  In  this 
case,  We  =  17.6,  and  e  =  0.001  for  both  simulations  presented  in  Fig.  10.  The  evolution  of  a  jet  under  steady 
back  pressure  (a)  is  compared  with  a  simulation  in  which  the  orifice  pressure  drop  is  varied  sinusoidally  about 
the  mean  (eqn  22)  at  a  frequency  of  Ar  =  2  (b).  The  unsteady  response  is  not  dramatically  different  than 
the  steady  case  even  under  the  large  amplitude  (75%)  perturbation  in  gas  pressure.  In  this  case,  the  orifice 
channel  provides  substantial  attenuation  of  the  imposed  oscillation  in  the  downstream  pressure. 


Conclusions 

A  stable,  time-accurate,  integration  scheme  has  been  developed  for  the  coupled,  nonlinear  dynamics  of  a 
free  surface  in  a  two-fluid  flow.  The  proceedure  is  illustrated  through  a  series  of  axisymmetric  and  2-D  flow 
examples  involving  liquid  jets  subjected  to  significant  pressure  disturbances  from  a  surrounding  gas.  While 
the  accuracy  of  the  time  integration  becomes  more  critical  in  cases  where  stagnation  points  are  present  on 
the  free  surface  the  technique  is  shown  to  work  well  even  in  flows  where  stagnation  points  are  present.  As 
the  gas/liquid  density  ratio  approaches  unity,  the  stability  and  accuracy  of  the  schemes  degrade.  However, 
the  methods  described  herein  have  been  applied  to  physically-meaningful  gas/liquid  flows  which  generally 
exhibit  density  ratios  less  than  0.01  even  under  high  pressure  conditions. 
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Nomenclature 

a  =  orifice  or  column  radius 
G  =  free  space  Greens  function 
k  =  wave  number 

K{jf>)  =  complete  elliptic  integral  of  the  first  kind 
/io.i  =  modified  Bessel  functions  (eqns  16  and  17) 
n  —  coordinate  normal  to  local  surface 
P  =  pressure 

q  =  velocity  normal  to  local  boundary 
r  =  radial  coordinate  (axisymmetric  flow) 
s  =  coordinate  aligned  with  local  surface 
t  =  time 

We  =  Weber  number,  We  —  pgU’^ajcr 
X  =  axial  coordinate  (2-D  flow) 
y  =  transverse  coordinate  (2-D  flow) 

2  =  axial  coordinate  (axisymmetric  flow) 
a  “  boundary  point  singularity  contribution 
(3  =  surface  slope 
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c  =  gas/liquid  density  ratio 
r  =  domain  boundary 
K  =  surface  curvature 
0  =  velocity  potential 
p  =  density 
cr  =:  surface  tension 
u;  =  frequency 

Subscript 
(  )^  =  gas  phase 
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Figure  1.  Schematic  of  Computational  Domain  for  Liquid  Jet  in  Wind-Induced  Regime  Denoting  Bound¬ 
ary  Conditions 

Figure  2.  Nonlinear  Jet  Evolution  in  First  Wind-Induced  Regime 
Figure  3.  Nonlinear  Jet  Evolution  in  the  Second  Wind-Induced  Regime 

Figure  4.  Computational  Domain  and  Boundary  Conditions  for  Liquid  Column  Subjected  to  Acoustic 
Perturbations 

Figure  5.  Time  History  for  Top  Node  Position;  We  =  0.1,  c  =  0.01, 

Figure  6.  Column  Shape  at  Various  Times  During  Excitation  Process 

Figure  7.  Density  Ratio  Influence  on  Maximum  Aspect  Ratio  of  the  Column 

Figure  8.  Schematic  of  Computational  Domain  for  Liquid  Jet  Injected  into  Quiescent  Gas 

Figure  9.  Effect  of  Gas  Density  on  Initial  Liquid  Jet  Behavior,  We  —  17.6 

Figure  10.  Effect  of  Unsteady  Chamber  Conditions  on  Jet  Evolution,  We  —  17.6,  e  =  0.001,  (a)  -  Steady 
Flow;  (b)  -  75%  Pressure  Oscillation  About  Mean  Flow  a,i  k  =  2 
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Figure  1:  Schematic  of  Computational  Domain  for  Liquid  Jet  in  Wind-Induced  Regime  Denoting  Boundary 
Conditions 


Figure  2:  Nonlinear  Jet  Evolution  in  First  Wind-Induced  Regime 


Figure  3:  Nonlinear  Jet  Evolution  in  the  Second  Wind-Induced  Regime 
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Figure  4:  Computational  Domain  and  Boundary  Conditions  for  Liquid  Column  Subjected  to  Acoustic  Per 
turbations 
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Figure  7:  Density  Ratio  Influence  on  Maximum  Aspect  Ratio  of  the  Column 


Jet  Orifice 

Figure  8:  Schematic  of  Computational  Domain  for  Liquid  Jet  Injected  into  Quiescent  Gas 
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t  =  8 


(a)  T  Const.  AP 

(b)  -  Osc.  AP 


Figure  10:  Effect  of  Unsteady  Chamber  Conditions  on  Jet  Evolution,  We  =  17.6,  (a 
^  75%  Pressure  Oscillation  About  Mean  Flow  at  k  =  2 
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)  -  Steady  Flow;  (b)  - 
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